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ABSTRACT 

The Analytic Predictor Corrector (APC) and Energy Controller (EC) atmospheric 
guidance concepts have been adapted to control an interplanetary vehicle aerobraking in 
the Martian atmosphere. Modifications are made to the APC to improve its robustness to 
density variations. These modifications include adaptation of a new exit phase algorithm, 
an adaptive transition velocity to initiate the exit phase, refinement of the reference dy- 
namic pressure calculation and two improved density estimation techniques. The modi- 
fied controller with the hybrid density estimation technique is called the Mars Hybrid 
Predictor Corrector (MHPC), while the modified controller with a polynomial density esti- 
mator is called the Mars Predictor Corrector (MPC). 

A Lyapunov Steepest Descent Controller (LSDC) is adapted to control the vehicle. 
The LSDC lacked robustness, so a Lyapunov tracking exit phase algorithm is developed to 
guide the vehicle along a reference trajectory. The equilibrium glide enti^ phase is em- 
ployed for the first part of the trajectory. This algorithm, when using the hybrid density es- 
timation technique to define the reference path, is called the Lyapunov Hybrid Tracking 
Controller (LHTC). With the polynomial density estimator used to define the reference 
trajectory, the algorithm is called the Lyapunov Tracking Controller (LTC). 

These four new controllers are tested using a six degree of freedom computer simu- 
lation to evaluate their robustness. MARS-GRAM is used to develop realistic atmo- 
spheres for the study. The atmospheres are then perturbed using square wave density 
pulses. The MHPC, MPC, LHTC and LTC show dramatic improvements in robustness 
over the APC and EC. The MHPC, MPC, LHTC and LTC all complete the initial phase of 
testing (using square wave density pulses) with no failures. The second phase tests the 
MHPC, MPC, LHTC and LTC against atmospheres where the inbound and outbound den- 
sity functions are different. Square wave density pulses are again used, but only for the 
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outbound leg of the trajectory. Additionally, sine waves, in both altitude and range, are 
used to perturb the density function. All four new controllers are able to compensate for 
the outbound leg density pulses with no hard failures, but the algorithms are sensitive to 
large amplitude density pulses. Additionally, these control algorithms are sensitive to 
large amplitude sine waves, particularly sine waves in range. The hybrid density estimator 
responds poorly to sine waves in range with wavelength between twenty and two hundred 
nautical miles. The polynomial density estimator is sensitive to wavelengths between five 
hundred and two thousand nautical miles. Overall, the polynomial density estimator per- 
forms better than the hybrid density estimator. The Lyapunov tracking phase performs 
better than the predictor correctors and the LTC is the most robust control algorithm exam- 


ined. 
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CHAPTER I 
INTRODUCTION 

When orbital transfer is required near a celestial body with an atmosphere of suffi- 
cient altitude and density, it is often advantageous to utilize aerodynamic forces to aid in 
the transfer 1 ' 11 . Aerodynamic drag forces are used to reduce the kinetic energy, while 
aerodynamic lift forces are used to control the trajectory during the maneuver. The result 
is a vehicle weight savings equivalent to the propellant necessary to perform the maneuver. 
The critical factor for success in the aerobraking maneuver is the performance of the guid- 
ance control system. The National Aeronautics and Space Administration plans a 1992 
launch of the Aeroassisted Flight Experiment (AFE) to serve as a proof-of-concept and 
test vehicle for aerobraking orbital maneuvers . Meanwhile, an aeroassisted orbital 
transfer maneuver is planned for the Mars Rover/Sample Return (MRSR) Mission to re- 
duce the orbit energy from the hyperbolic Martian approach orbit to capture into a low 
Mars orbit with a commensurate AV savings of over 8000 ft/s when compared with an all 
propulsive mission . 

Fig. 1 shows the proposed interplanetary mission concept. The vehicle will be 
launched from Earth into an elliptical heliocentric orbit. The vehicle will travel almost 
eight months in this interplanetary orbit making mid-course corrections as necessary to in- 
tercept Mars. When the vehicle reaches Mars there will be 6 km/sec difference between 
the velocity of the vehicle and Mars orbital velocity. Without some method of changing 
the vehicle’s velocity it would swing by Mars without capturing into a Martian orbit. The 
proposed method for imparting this velocity change is to use the aerodynamic forces im- 
parted by the Martian atmosphere. The final mid-course correction to the interplanetary 
orbit will allow the vehicle to enter the Martian atmosphere as shown in Fig. 2. The typi- 


Joumal model is AIAA Journal of Aircraft. 
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cal sequence of events for an aerobraking maneuver call for the vehicle to plunge into the 
atmosphere and fly deep in the atmosphere until the velocity is appropriately reduced. The 
vehicle then executes a pullout maneuver exiting the atmosphere in a low Mars orbit. Fi- 
nally, a series of propulsive maneuvers are performed to transition the vehicle into the de- 
sired final orbit. 



Fig. 2 Aerobraking Maneuver Sequence of Events 
(Adapted from Reference 14) 
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In the past, for space missions reentering the Earth’s atmosphere, only the destina- 
tion coordinates have been specified. In targeting to the correct orbit following the aero- 
braking maneuver the guidance system must accurately control the final position of the 
vehicle as well as the final velocity vector. The atmospheric lift and drag forces affecting 
the vehicle are proportional to the atmospheric density, but atmospheric density is highly 
variable 15 ' 22 . The guidance algorithm must be robust enough to control to the final state 
even with these uncertainties. The focus of this report is to study the relative merits of 
several existing and novel guidance algorithms, with particular emphasis upon the extent 
to which the algorithms tolerate our ignorance of the Martian atmosphere. 

Although robustness with respect to density variations was a prime factor in develop- 
ing and choosing the guidance scheme for the AFE 3 ' 10 ' 11,23 " 25 it becomes even more criti- 
cal for the Mars mission. Scientists have worked for many years to characterize the 
Earth’s atmosphere. Accelerometer data gathered during space shuttle returns have al- 
lowed us to characterize not only the average density values but also the expected magni- 
tude and frequency of the random density variations . In designing the AFE guidance 
system the Earth’s upper atmosphere in the region 250000-400000 feet, was assumed to 
have density variations of ±25% from standard values over small altitude intervals 15 . The 
Martian atmosphere goes through global atmospheric expansions and contractions equiva- 
lent to an atmospheric shift of 10 km 22 . Additionally, data gathered from the Viking I and 
Viking II landers show density variations of 20 to 30% over small altitude intervals in the 
aerobraking region 18 ' 19 . Since we have only two sets of density measurements from the 
Martian aerobraking region, we must expect even larger density variations to occur. It is 
conjectured that density shears of 60% or greater may be encountered. Development of a 
robust controller capable of acceptable performance given the large, unpredictable density 
variations in the Martian atmosphere is, therefore, vital to the success of the MRSR mis- 


sion. 
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Background 

The basic technology required to perform hypersonic flight in an atmosphere was de- 
veloped in the 1950s to support the development of intercontinental ballistic missiles 26 . 
Technology was extended to allow the Mercury and Gemini projects to dissipate kinetic 
energy by entering the atmosphere with low ballistic coefficient vehicles. Major advances 
in entry technology were made during the Apollo program, especially in the areas of navi- 
gation, guidance and control during atmospheric maneuvering. With the Space Transpor- 
tation System (STS) came a reusable capability to deploy and retrieve satellites from Low 
Earth Orbit (LEO). Deployment of the Space Station will provide a permanent base in 
LEO capable of performing maintenance and repair of satellites. However, with a large 
percentage of satellites in Geosynchronous Earth Orbit (GEO)an economical system of 
deploying satellites to GEO and then returning them to LEO is required. The National 
Aeronautics and Space Administration (NASA) has developed an aerobraking vehicle, the 
AFE, to meet the return requirement 12 . In designing the Mars Rover/Sample Return mis- 
sion an aerobraking phase similar to that of the AFE is envisioned to dissipate kinetic en- 
ergy from the hyperbolic Martian approach orbit leaving the satellite in a Low Mars 
Orbit 13 . 


The AFE, scheduled for launch in 1992 will serve as a proof of concept and test ve- 
hicle for aerobraking. AFE will enter the atmosphere with the same velocity as a vehicle 
returning from GEO. The vehicle will fly trimmed at a constant angle of attack, and there- 
fore, at near a constant lift to drag (L/D) ratio. AFE will roll about the velocity vector to 
modulate the in plane portion of the lift to control the trajectory while drag dissipates ki- 
netic energy. The vehicle will exit the atmosphere, after an appropriate energy reduction, 
with exit velocity and flight path angle such that the final orbit will rendezvous, at apogee, 
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with the desired LEO. The mission concept for the Mars aerobrake is expected to be quite 
similar, but of course, for the MRSR objectives. 

The complexity of the competing inequality and equality constraints placed on an 
aerobraking maneuver make definition of an optimal, robust control algorithm extremely 
difficult. The simplest controllers are open loop controllers designed to optimize the tra- 
jectory for a specific atmosphere, entry conditions and vehicle design. Talay, et al, opti- 
mized a bank angle history for a nominal 1962 atmosphere using a trajectory optimization 
code. When this bank angle history was used in trajectory simulations with off-nominal 
atmospheres in several cases the vehicle either exited the atmosphere early or failed to exit 
at all. Vinh 4 first formulated an optimal , minimum fuel, control problem using a com- 
bined propulsive and aerodynamic transfer. He shows that an optimal combined propul- 
sive/aerodynamic orbit transfer will require only 32% of the total AV required for an all 
propulsive maneuver for an orbit transfer from GEO to LEO. Then Vinh, et al, 27 produce 
an explicit guidance scheme for the aerobraking phase of a drag modulated aeroassisted 
transfer between elliptical orbits. They find the optimal strategy consists of bang-bang 
control but then point out that the strategy is difficult to realize because the switching time 
must be very accurate, “within a fraction of a second to avoid crashing.” They propose an 
alternative strategy whereby the drag is controlled between minimum and maximum val- 
ues as a function of the current state. Kechichian, et al, 28 also acknowledge that for a drag 
modulated vehicle bang-bang control is optimum for minimizing the total A V required to 
achieve the desired orbit, but in an effort to reduce the sensitivity to switch point timing a 
new controller is developed to add an additional degree of control. 

Sensitivity analysis shows that this control scheme has essentially zero sensitivity to an at- 
mospheric density profile of ±15% of nominal but an entry corridor width of ±0.1° 
should be maintained to avoid excessive A V requirements. 
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Much work has been performed in the area of optimal aeroassisted plane changes. 
Hull, et al, 29 derives an optimal guidance scheme for performing an aeroassisted plane 
change between circular orbits. They assume a parabolic drag polar for the vehicle and 
use Loh’s constant 30 to include gravitational terms and apparent lift terms in the analysis. 
They find bank angle and angle of attack time histories which minimize the total AV re- 
quired to perform the maneuver by maximizing the exit velocity following the aero phase. 

o 1 

Plane changes of 10 to 40 degrees are demonstrated. Later the problem is reformulated 0 
using heading as the independent variable and assuming that Loh’s term may be either 
positive or negative. They show that only one solution exists and that it may be found by 
solving a fourth order polynomial. Hull, McClendon, and Speyer 32 then reform the prob- 
lem assuming an elliptic drag polar and obtain similar results. They show that near the 
end of the atmospheric turn Loh’s term is not constant which may cause extremely high 
angles of attack. Finally, Hull and co-workers 33 assume Loh’s term is piecewise constant 
during the turn and reformulate the problem. Using the method of successive approxima- 
tions they construct a control law which results in a final velocity within 1 % of the true op- 
timal final velocity for a 40° plane change and results in a very reasonable maximum 
angle of attack of 30°. Johannesen, et al, 5 formulate an approximate control law for lift 
and bank angle to maximize orbit plane change using an aeroassist maneuver. The control 
law is tested for a wide range of speed ratios V g /Vj. They observe that the maximum turn 
angle for any speed ratio is proportional to the maximum lift-to-drag ratio. 

Two unique methods of determining atmospheric guidance control laws have been 
developed. Mease and McCreary 6 propose using an approximate closed form solution of 
the equations of motion. Their solution divides the trajectory into three regions. During 
the beginning and end of the trajectory the gravitational terms are assumed to dominate, 
while in the mid-portion of the trajectory aerodynamic terms are assumed to dominate. 
The solutions for each of these regions is combined using the method of matched asymp- 
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totic expansion. Final apocenter values within 9% of the targeted values are demonstrated 
for a wide range of entry flight path angles for a simulated Mars aerocapture mission. The 
other unique method developed by Lee and Grantham^ uses Lyapunov optimal feedback 
control to minimize the A V required to raise perigee following an aerobraking maneuver. 
This method calculates a descent function and then seeks to move the system in a pre- 
ferred direction, opposite the gradient of the descent function. The Lyapunov feedback 
controller is compared with an optimal open loop controller derived using calculus of vari- 
ations for the nominal 1962 standard atmosphere. Superior, robust performance for the 
Lyapunov controller is demonstrated for both the standard atmosphere and a shuttle-de- 
rived atmosphere. 

Control laws developed using optimal control theory offer excellent performance in 
numerical simulations, but those methods which require extensive computation for each 
control update have been at a distinct disadvantage due to limitations of onboard comput- 
ing capability. For this reason several simplified guidance schemes have been developed. 

o 

Letts and Pelekanos developed a control law using bank -angle modulation of the lift vec- 
tor to establish a constant axial deceleration level until the required exit velocity is 
reached, when full lift up is commanded. They show that AV required to circularize fol- 
lowing the aerobrake maneuver increases approximately 35 m/s for each percent change 
from the nominal value for an atmosphere that is multiplied by a constant density bias. 
Gamble, et al, develop a control scheme similar to that of Letts and Pelakanos except 
Gamble’s method commands to an equilibrium glide rather than a constant axial decelera- 
tion until the desired velocity is achieved. After the desired velocity is achieved, full lift 
up is again commanded for atmospheric exit. Gamble finds that a 50% increase in density 
has little effect on the total AV required but a 50% decrease in density increased A V re- 
quired by about 35% due to problems in establishing the equilibrium glide. Cerimele, et 

A 

al, use Gamble’s equilibrium glide phase during the entry portion of the trajectory, but 
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then switch to a reference drag profile like Letts and Pelekanos for the exit portion of the 
trajectory. Density shears in the atmosphere are simulated and the AV required following 
the aerobraking maneuver is found to be very sensitive to density ratios exceeding ±30% 
occurring over altitude ranges of 1,000 to 10,000 ft. Cerimele and Gamble 9 produce an 
analytic predictor corrector guidance algorithm, again using the equilibrium glide entry 
phase but with a predictor corrector exit phase designed to target apogee more accurately. 
The predictor corrector algorithm assumes a constant altitude rate and an exponential at- 
mosphere to predict apogee. The predictor algorithm iterates altitude rate until a value is 
found which produces the desired apogee. The vehicle is then commanded to this altitude 
rate. An interesting feature added to this algorithm is a low pass density filter. Density is 
computed on-board based on accelerometer data. Calculated density is then compared 
against predicted density values and future predicted values are adjusted accordingly. This 
guidance algorithm was tested numerically using combined dispersions of ±0.2° in entry 
flight path angle, ±20% density variations and ±33% L/D. The final apogee value was 
within 2 nm of the target value in all cases. Gamble, et al, 34 present three atmospheric 
guidance concepts for aeroassist orbit transfer vehicles. The first method presented is the 
Analytic Predictor Corrector, already discussed. The second is a Numerical Predictor 
Corrector algorithm which numerically integrates a trajectory assuming constant bank an- 
gle magnitude and an assumed density profile. The bank angle is iterated until the desired 
apogee is computed and the vehicle is commanded to this bank angle. The final control al- 
gorithm presented is the Energy Controller which guides the vehicle to a desired energy 
state at atmospheric exit. The energy gain, defined as the ratio of energy rate to energy er- 
ror, is controlled so that energy error exponentially goes to zero at atmospheric exit. The 
energy gain command is converted to an altitude rate command which in turn is converted 
to a bank angle command. Their results show that all three algorithms are capable of 
maintaining the final apogee within 10 nm and AV within 50 ft/sec of the nominal values 
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for test cases with dispersions of ±4 nm in perigee, ±50% in density, ±50% in W/CqA 
and +50% in L/D. The analytic predictor corrector and Energy Controller show slightly 
worse results for the -50% L/D case. 

Fitzgerald and Ward 10 ' 11 investigate the sensitivity to density shears of the Analytic 
Predictor Corrector and Energy Controller algorithms while guiding the AFE vehicle. 
They consider spike and step shaped density dispersions of ±10 and ±20% magnitude 
with durations of 5,000 and 10,000 feet, starting at altitudes between 260K and 295K ft. 
A V increases up to 60% for the Energy Controller and 41% for the APC are demonstrated. 
Fitzgerald 1 1 then formulates a Hybrid Predictor Corrector algorithm which uses the atmo- 
spheric density profile determined during the entry phase in the predictor corrector of the 
exit phase. This significantly reduces the sensitivity to density shears for atmospheres 
where the exit atmosphere matches the entry profile. 

Meyerson and Cerimele review the aeroassist vehicle requirements for the Mars 
Rover/Sample Return Mission. They use a modified analytic predictor corrector algorithm 
referred to as HYPAS as the controller for vehicles with L/D ranging from 0.3 to 1 .5. Ad- 
ditionally, entry velocities from 5.79 to 9.20 km/sec were investigated. A recommenda- 
tion of this study is, “to refine the HYPAS guidance algorithm to control the trajectory 
more accurately in the exit phase.” They recommend using two exponential atmosphere 
models in the guidance predictor. 

“The Mars atmosphere is highly variable on a daily, seasonal and annual basis 18 .” 
The thin atmosphere and solar heating produce a large daily temperature range which 
translates to a large daily density fluctuation 18 ' 21 . Fig. 3 shows that at the surface the 
Martian atmospheric density is approximately two orders of magnitude less than that of 
the Earth’s atmosphere and at aerobraking altitudes there is still more than an order of 
magnitude difference between the density of Earth and that of Mars. “On an annual basis, 
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Fig. 3 Comparison of COSPAR Northern Hemisphere Mean Mars 
Density Model and 1962 U.S. Standard Earth Atmosphere 
(Adapted from Reference 19) 


the atmospheric pressure at the surface changes by ±15% due to condensation and subli- 
mation of the CO 2 17 ” which produces a global expansion and contraction of the atmo- 
sphere of roughly 10km. Fig. 4 presents the density deviations of the COSPAR high 
density model and the COSPAR low density model relative to the COSPAR Northern 
Hemisphere Mean Model. Global dust storms absorb radiation high in the atmosphere, 
thereby increasing the upper atmosphere temperature and causing a large scale expansion 
of the atmosphere. The density is then substantially increased at orbital and entry alti- 
tudes. Additionally, density of the Martian atmosphere varies widely on a daily basis. 
Fig. 5 shows the expected morning and afternoon density profiles calculated for summer at 
the Viking 1 lander location while Fig. 6 shows the calculated density profiles for winter at 
the Viking 1 lander location. These figures show that at aerobraking altitudes the density 
may vary by as much as 100 to 150% on a daily basis. The Viking 1 and 2 landers mea- 
sured atmospheric properties during their descent and recorded peak to peak density varia- 
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Fig. 4 Density Deviations of the COSPAR Low-Cool and COSPAR 
High-Warm Density Models as Compared to the COSPAR Northern 

Hemisphere Mean Model 
(Adapted from Reference 19) 


tions in the aerobraking region of 30% over a 15 km altitude band and 20% over a 10 km 
region ‘ . Fig. 7 and Fig. 8 present the density variations measured by Viking 1 and Vi- 
king 2 respectively during their descent to Mars. The Mars Global Reference Atmosphere 
Model (MARS-GRAM) 35 characterizes Mars atmospheric properties, density, tempera- 
ture, pressure and wind speed and direction, as functions of date, time, latitude, longitude, 
altitude, solar activity and dust storm activity. 

This report will characterize the sensitivity of selected aerobraking guidance algo- 
rithms with respect to density variations of the type and magnitude expected in the Mar- 
tian atmosphere to determine their suitability to perform the MRSR Mission. A guidance 
algorithm capable of acceptable performance in spite of the uncertainty in Martian atmo- 
spheric density or methods of reducing the uncertainty will be developed. To attain this 
goal the following research objectives are proposed. 
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Fig. 5 Morning and Afternoon Density Profiles Calculated for the 
Viking 1 Lander Location During the Summer 
(Adapted from Reference 19) 



Fig. 6 Morning and Afternoon Density Profiles Calculated for the 
Viking 1 Lander Location During the Winter 
(Adapted from Reference 19) 
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Research Objectives 

When designing a control law for an aeroassist maneuver, an exponential variation 
of atmosphere density with altitude is an extremely attractive computational simplifica- 
tion. However, given the large density biases and density shears of the Martian atmo- 

10 

sphere ' , a guidance algorithm optimized for the MRSR vehicle aerobraking in the 
assumed exponential atmosphere may demonstrate poor performance and potentially cata- 
strophic failures if realistically off-nominal conditions are encountered 13 . Especially 
when errors in navigation accuracy and/or vehicle L/D are also considered, the results of 
using any fixed density model may be catastrophic, with the vehicle either not being cap- 
tured into a Martian orbit or failing to exit the Martian atmosphere. As a result, the sensi- 
tivity of MRSR atmospheric guidance to perturbations in density as well as to navigation 
errors and L/D errors is critical to the success of the mission. 

A systematic method of evaluating the effects of density biases and density shears in 
combination with navigation errors and L/D errors on an MRSR atmospheric guidance al- 
gorithm is sought. The methods established will be used to evaluate candidate guidance 
algorithms, including algorithms developed in this report. Toward these objectives, the 
first task is to determine the expected extremes in Martian density. MARS-GRAM 35 will 
be utilized for this task. The highest and lowest density atmospheres expected are deter- 
mined as a function of season, time of day, solar activity and dust storm activity. Since 
there have only been two space probes which have measured density through the Martian 
aerobraking region, the nature and magnitude of expected worst case density shears is not 
known precisely and must be estimated. These atmosphere extremes are checked against 
the Mars standard atmospheres 18 and the Viking data 18 ' 19 . There is a proposal for the 
Mars Aeronomy Observer (MAO ) 20 ' 21 to send additional probes to the Martian surface to 
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gather more data quantifying these shears; but, this mission is still years away. The ex- 
pected navigational accuracy and probable L/D errors are also be determined. 

Secondly, the Analytic Predictor Corrector algorithm, 9 the algorithm chosen for the 
AFE mission, and the Energy Controller 34 are fine tuned for the MRSR mission. Then the 
sensitivity of these algorithms when faced with these density biases, density shears, navi- 
gation and L/D errors are determined. The six degree-of-freedom Program to Optimize 
Simulated Trajectories (POST) 36 is used in this analysis. The sensitivity of these algo- 
rithms is visually presented by plotting three dimensional sensitivity surfaces with the 
qualitative objective of finding the worst combinations of dispersions and defining the per- 
formance bounds of these two controllers. Methods of improving the performance of 
these algorithms, especially methods of using information derived early in the trajectory, 
to improve the performance in the latter portions of the trajectory (similar to the methods 
proposed by Fitzgerald in his Hybrid Predictor Corrector algorithm) are evaluated. Two 
new algorithms called the Mars Hybrid Predictor Corrector (MHPC) and the Mars Predic- 
tor Corrector (MPC) are developed. These two algorithms differ only in their density esti- 
mation techniques. This task, along with developing the new algorithm proposed below, 
are crucial to the research and secondary only to the task of defining absolute robustness 
limits. 

The third order of business is to explore more elegant ways to optimize the control- 
ler, especially ways of improving the robustness. A potential candidate Lyapunov Steep- 
est Descent Controller 7 (LSDC) similar to the one suggested by Lee and Grantham is 
coded and its performance tested against the same perturbations as the others. Two new 
algorithms are developed, again differing only in their density estimation techniques. 
They are called the Lyapunov Tracking Controller (LTC) and the Lyapunov Hybrid Track- 
ing Controller (LHTC). 
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Finally, and most important of all, the robustness limits of the improved MPC, 
MHPC, LTC and LHTC controllers are characterized. The guidance performance is thor- 
oughly tested to find the tolerable limits on density bias and density shear given the proba- 
ble errors in navigation and L/D. POST is used to test the guidance algorithms, using the 
Viking atmosphere profiles 18 ' 19 . The limits on atmosphere dispersions, considering the 
inherent navigation and probable LTD errors, under which acceptable controller perfor- 
mance will occur is clearly defined from the results of these simulations. 

These limits are checked against the worst case perturbations expected for the mis- 
sion ' * . Conclusions are drawn regarding the performance of these algorithms when 
faced with the expected density variations, as well as possible variations in vehicle lift to 
drag ratio and entry flight path angle. Recommendations for future study are then presented. 

Organization of the Report 

Improvements made to the Analytic Predictor Corrector and Hybrid Predictor Cor- 
rector control algorithms are presented in Chapter II. Derivation of the LSDC and a LTA 
are presented in Chapter III. Chapter IV details the model used for the trajectory simula- 
tions and the atmospheric models. Vehicle mass and aerodynamic data are presented 
along with atmospheric entry conditions. The controller test program is outlined and the 
perturbations in atmospheric density, vehicle lift and drag, and entry conditions which 
were used in the test program are presented. Finally the performance of the various con- 
trollers is presented. In Chapter V the four best performing controllers, the MPC, MHPC, 
LTA and the LHTA algorithms, are tested against each other in a head to head fashion. 
The perturbation limits which define the edge of the envelope where acceptable perfor- 
mance is attainable are determined. Conclusions and recommendations are presented in 
Chapter VI. 
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CHAPTER II 

IMPROVEMENTS TO THE 
PREDICTOR CORRECTOR ALGORITHM 

The Analytic Predictor Corrector controller, discussed by Cerimele and Gamble 9 , is 
the control algorithm selected for the AFE 12 . This controller was adapted to the MRSR 
problem and tested against expected perturbations in the Martian atmosphere as well as 
perturbations in entry conditions and vehicle lift and drag characteristics. While testing 
the Analytic Predictor Corrector controller, several areas were found where improvements 
could be made. The constant multiplier used to determine the reference dynamic pressure 
was changed in an effort to gain robustness and prevent premature exit from the atmo- 
sphere. Borrowing a technique first employed by Fitzgerald 11 , an improved atmospheric 
model used by the predictor step to determine velocity loss during the exit phase was also 
incorporated. Then a new method of estimating density incorporating a polynomial to fit 
the normalized density function was developed 42 . A modified exit phase, first developed 
at the Charles Stark Draper Laboratory 43,44 , was incorporated and tested, first without and 
then with the improved atmospheric models. The new exit phase also assumes a constant 
altitude rate to compute the velocity loss due to aerodynamic drag. However, instead of 
predicting the exit state and iterating altitude rate to target the desired apocenter, the veloc- 
ity required to attain the desired apocenter altitude is computed based on the current state 
and an estimate of the remaining velocity loss due to aerodynamic drag. The iteration is 
simplified to a single step altitude rate correction. With the large density shifts present in 
the Martian atmosphere and the uncertainties in vehicle and entry conditions the velocity 
at which the controller transitions from the equilibrium glide phase to the exit phase (in- 
corporated as a controller constant for the APC controller operating in the Earth atmo- 
sphere) was changed to an adaptive parameter. 
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The predictor corrector algorithm developed here which uses a variation of Fitzger- 
ald’s density estimation scheme is referred to as the Mars Hybrid Predictor Corrector 
(MHPC). The algorithm which incorporates the polynomial density estimator is called 
simply the Mars Predictor Corrector (MPC) algorithm. The modifications presented here 
convert the predictor corrector algorithm from a good controller for guiding the aerobrak- 
ing phase of a space vehicle returning from Geosynchronous Earth Orbit (GEO) into a ro- 
bust control algorithm capable of accurately guiding an interplanetary probe through an 
aerobraking maneuver in the Martian atmosphere. 


Reference Dynamic Pressure Calculation 


The equilibrium glide phase of the APC controller seeks an equilibrium condition 
with the vehicle following a reference dynamic pressure path. The reference dynamic 
pressure is calculated as a multiple, K-, of the dynamic pressure required to maintain 
equilibrium with the lift vector oriented down. Gamble, et al^, recommend that this mul- 
tiplier be 1 .33 for the AFE which aerobrakes in the Earth’s atmosphere. 
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Ideally, to have the minimum AV required after the aerobraking maneuver the veloc- 
ity of the vehicle should be decreased as high in the atmosphere as possible. Some studies 
have considered, as a minimum AV aeromaneuver, the case of a vehicle with infinite lift 
skimming the edge of the atmosphere until the velocity has decreased by the appropriate 
amount so the vehicle can be released into a Hohman Transfer orbit from the circular orbit 
at the edge of the atmosphere to the desired orbit^. However, decreasing the velocity high 
in the atmosphere means flying in a region of lower density and consequently lower dy- 
namic pressure. Flying higher requires the in-plane component of the lift vector to be ori- 
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ented more downward to maintain equilibrium. This idealization is satisfactory only with 
a smooth exponential density profile; density shears are not allowed. If the vehicle is fly- 
ing in equilibrium using 50% of the lift down capability (O = ±120°) and a density shear 
is encountered which decreases the density by 50% suddenly, full lift down will be re- 
quired to maintain equilibrium. In actuality, because of time lags between the encounter 
of a density shear and the vehicle’s response, coupled with the potential necessity of re- 
ducing a positive altitude rate, the minimum acceptable reference dynamic pressure to 
maintain control, if the density suddenly decreases by 50%, is considerably more than 
twice that required to maintain equilibrium in a full lift down configuration. 

One potential drawback to increasing K- is that the trajectory loads are increased 
over a portion of the flight. Heat rates and vehicle acceleration loads are increased for the 
portion of the portion of the trajectory after the minimum altitude point until the transition 
to the exit phase, however, for the range of K- between 1.33 and 4.5 the maximum heat 
rates, g loading, the minimum altitude, and even the maximum dynamic pressure do not 
change. 4.5 was the largest value which would not adversely affect the peak trajectory 
loads. Furthermore, it has been found that the total heat integrated heat load calculated is 
lower for a higher value of K- because the vehicle’s deceleration is greater and less time is 
required to reduce the vehicle’s velocity. 

Fig. 9 presents the altitude time histories and Fig. 10 presents the velocity time histo- 
ries for trajectories flown through a nominal Martian atmosphere perturbed with a square 
wave pulse of 20,000ft duration located between 140,000ft and 160,000ft altitude. This 
pulse multiplies the nominal density function by 0.5 in this altitude region. The three 
curves presented in each figure represent K- values of 1.33, 2.2 and 4.5. Notice that the 
trajectories where K- equals 1.33 and 4.5 perform well. The final apocenter is within 
three nautical miles of the targeted 270 nm for both cases. However, the trajectory flown 
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with K- equal to 2.2 skips out, or exits the atmosphere early while the velocity is still too 
high. Skip outs are difficult to predict because they are caused by a sudden negative densi- 
ty shear which reduces the vehicle’s available lift and drag force. If this negative density 
shear occurs when the vehicle is in a relatively safe regime where the sudden decrease in 
density will not place the vehicle in a critical situation, a skip out may not occur. If the ve- 
hicle has a positive altitude rate while reversing the bank with the lift vector oriented up, 
or if the control scheme allows the vehicle to overshoot the reference dynamic pressure 
trajectory, thus temporarily flying at a dynamic pressure lower than the reference value, a 
skip out is quite likely. Combinations of these factors are an even greater challenge for 
the controller when a sudden negative density shear is encountered. Increasing K- will 
not always prevent a skip out, but increasing K- does tend to reduce the probability of a 
skip out. Indeed, with K- set to 4.5 (the largest value possible without adversely affecting 
peak trajectory loads) no skip outs were encountered during the validing simulations. 


With the uncertainties in the Martian atmosphere the slight penalty in AV required to 
increase this multiplier to 4.5 (less than 10 ft/sec), when compared to the 1 .33 value rec- 
ommended for the Earth atmosphere, seems to be a small penalty to gain additional ro- 
bustness and limit the possibility of a premature skip out. q rtf is therefore calculated 
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Improved Exit Phase Density Models 

The second area of improvement is the density estimation technique. Good density 
estimation is critical for the success of any Martian aerobraking. With the wide density 
variations possible in the Martian atmosphere, the correct path, given an estimate of future 
density, may prove disastrous if that estimate is wrong. Given that the MRSR vehicle will 
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traverse over 1000 nautical miles during the aerobraking maneuver it is entirely conceiv- 
able that the density function encountered may vary as much as the 100 to 150% variation 
in density between morning and afternoon presented in Fig. 5 and Fig. 6. The density es- 
timation technique must not only develop a profile of density versus altitude, but must 
continue to update this estimate throughout the trajectory based on the latest accelerome- 
ter-based density measurement. Two methods for performing this task are suggested here. 

Hybrid Density Estimator 

One possible approach is to model the atmosphere like Fitzgerald did 1 1 . During the 
entry phase accelerometer-derived density is recorded near each 1000 ft altitude interval 
along with the altitude for each density measurement. Then, during the exit phase, the 
predictor step uses these measurements to predict the velocity loss due to atmospheric 
drag. One difference between our approach and that of Fitzgerald is the inclusion of a 
density multiplier derived from accelerometer-generated density measurements which 
continue to update the density estimated throughout the trajectory. 

The accelerometer-generated density measurements taken during the entry phase of 
the aerobraking maneuver are, quite likely, the best estimate of the atmospheric density 
function available for the exit phase of the trajectory. These measurements will be close, 
in both space and time, to the density for the remainder of the flight. They will indicate 
the general state of the atmosphere, that is whether the C0 2 is in a condensed or sublimat- 
ed state, and they will provide an indication of the vertical wave structure of the atmo- 
sphere. They will not provide any information on horizontal waves which may affect 
density on the outbound leg. To compensate for this latter shortcoming, a density multipli- 
er and a low pass filter like the one presented for the APC^’ -* 4 


were used. However, in- 
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stead of dividing derived density by the density predicted by an exponential function, the 
divisor is the density predicted from the measurements taken during the entry phase. 


-(h-h^/hS 

P model — ^ \ e 


( 3 ) 


where p j is the density which was measured at the lower edge of the current altitude band, 
h j is the altitude at which this measurement was taken, and hS is the scale height for the 
atmosphere computed between this density measurement and the next measurement ap- 
proximately 1000 ft higher. 
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The density multiplier is then computed by dividing the accelerometer derived densi- 
ty by the density predicted for the current altitude, using the density model derived during 
entry. The result is filtered using a low pass filter to remove high frequency density devia- 
tions which would have minimal effect on the post-aerobraking apocenter. 


*p = (l-/0* p + K(p/p model ) 
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To use this modified atmosphere in the predictor step, rewrite the expression relating 
rate of change of velocity and rate of change of altitude 48 : 
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This equation may be integrated assuming a constant altitude rate to give the velocity 
loss due to atmospheric drag between two arbitrary altitudes /tj and h 2 . 
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and, with the scale height as previously calculated 
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This equation gives the relative velocity at h 2 as a function of the relative velocity at 
/i j and the densities and altitudes at the two locations. The method for employing this fea- 
ture in the predictor step of the control algorithm is to first use the velocity, density and al- 
titude at the current satellite location as the subscript 1 variables and to predict the velocity 
at the next interval where density measurements were stored during the entry using that al- 
titude and that density multiplied by the density multiplier discussed earlier as the sub- 
script 2 variables. Then that velocity may be used to compute the velocity at the next 
altitude band using the lower stored density and altitude values as subscript 1 variables and 
the next higher density and altitude measurements as subscript 2 variables. Notice that the 
density multiplier, when multiplied by each of the stored density measurements, will can- 
cel in all but one location. 



CK pPj ~ il 2^ (p 2 'j 
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This procedure is repeated until the exit relative velocity is computed. The velocity 
change expected between the current location and atmospheric exit may be calculated by 
subtracting the current relative velocity from the predicted exit relative velocity. 

AV= y rx -V r (10) 


Polynomial Density Estimator 

The second method of density estimation curve fits a sixth order polynomial in alti- 
tude to the normalized density function. This technique uses accelerometer derived densi- 
ty measurements at three trajectory locations to define a two phase exponential function. 
Derived density is recorded at one second intervals and then normalized by the exponen- 
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tial function in an effort to remove the underlying predominant exponential component. 
As the satellite reaches the bottom of the trajectory a batch estimate 42 is used to perform a 
weighted least squares fit of the polynomial coefficients to the resulting normalized func- 
tion. After that, a sequential estimate 42 is used to continue updating the coefficients of the 
polynomial for the remainder of the trajectory. 

Based on MARS-GRAM^ 5 generated data a two phase exponential function was 
chosen to normalize the density data. The underlying exponential component is assumed 
to be two exponential functions divided at 250,000 ft altitude such that the normalizing ex- 
pression p is 


-(h- 250000) /hS2 

p } e 

(h> 250, 000ft) 

-(h- 250000) /hSl 

p } e 

(/t < 250, 000ft) 


hSl and hS2 are the scale heights below and above 250,000 ft. p ; is the density at 
250,000 ft, determined using accelerometer derived density which is filtered using a low 
pass filter like the one presented in Eq. (5). The scale height hSl is found by using the fil- 
tered density measurement when the vehicle’s altitude rate first becomes positive and that 
at 250,000 ft in Eq. (4). similarly, hS2 is found using the measured density at 400,000 ft 
and the measured filtered value from 250,000 ft. The density value chosen at 400,000 ft is 
not the filtered version because at this early point in the trajectory the density filter has not 
had sufficient data to converge to a reliable estimate. 

After the altitude rate first becomes positive and the constants of Eq. ( 1 1 ) have been 
determined, the density values which were saved at one second intervals during the de- 
scent into the atmosphere may be normalized. The resulting data is fit with a sixth order 
polynomial in normalized altitude using a weighted least squares (batch) criterion to select 
the coefficients for the polynomial. A ninth order polynomial was originally chosen be- 



27 


cause the Viking 1 and Viking 2 atmospheric descent data (Fig. 7 and Fig. 8) shows six 
and five major density extremes respectively in the aerobraking region and the Viking 1 
data shows four additional local extremes. These data suggest that at least a seventh order 
polynomial is required to model even the major extremes accurately. Because computa- 
tional requirements for the density filter increase approximately as the square of the order 
of the polynomial, it was desired to use as low order polynomial as practical. After several 
iterations a ninth order polynomial appeared to be numerically ill-conditioned. A sixth or- 
der polynomial behaved much more consistently and was adequate for modeling the ex- 
pected density function. The density function is approximated by 

p (h) = p (h) [c j + c 2 x + c 3 x 2 + c 4 x 3 + c^x 4 + c 6 x 5 + c 7 x 6 ] (12) 

h 

where x is normalized altitude, x = -r- . 

h e 

The coefficients c } through c 7 are initially determined by a weighted least squares 
batch estimate presented by Junkins 42 . The procedure is to begin with the batch of m nor- 
malized density measurements 


yj = 



(13) 


taken at the m known altitude locations {h-) at m one second intervals until the altitude 
rate becomes positive. The altitudes are normalized by the atmospheric interface altitude 
to determine the x-s. The batch estimator must select the coefficients c l through c 7 so 
that 


yj = C J + c 2 Xj + c 3 Xj + c 4 Xj +c 5 Xj+ c 6 Xj + c 7 x 6 - +ej (1 4) 

where ej is the residual errors after selection of the coefficients. This equation may be 
written in matrix form 


Y = At+E 


(15) 



28 


where £ is the estimate of the coefficients Cj through c ? and A is the matrix 
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The batch estimator is tasked with selecting £ to minimize the weighted quadratic 
function of residual errors 


<p = e t we 


(17) 


where the weighting matrix selected is a diagonal matrix of weights applied to the residual 
error of each measurement. 


W = 


w } 0 0 

0 w 2 0 


lo 0 wj 


(18) 


Because the vehicle is traveling into a region of higher density which has greater 
impact on the satellite trajectory than does the thinner atmosphere near entry and exit and 
because more recent data were deemed to be more representative of future density than 
was older data, the weights were chosen to increase with time. An exponentially increas- 
ing weighting function was chosen which doubled the weighting after 1000 seconds. This 
weighting function was selected through experimentation which showed that a slower in- 
creasing weighting function did not respond quickly enough to abrupt density shears to 
produce adequate controller performance, while weighting functions that increased faster 
tended to ignore data gathered early in the trajectory and produced a poor estimate of den- 
sity in the upper altitude regions. This poor estimate had a distinctly negative impact on 
controller performance. The weights selected were 
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Equation (15) is solved for E and the result substituted into Eq. (17). After some ma- 
nipulation $ is expressed 

<p = Y T WY-2Y T WA£+£ T A T WA£ (20) 

To minimize <J> it is necessary that 

V^cp = -2A T WY + 2A T WAC = 0. (21) 

This equation is solved to obtain the weighted least squares normal equation for £ 

£ = (A 1 WA)~ X A t WY. (22) 

After the satellite altitude rate becomes positive the estimator switches to a linear se- 
quential estimator 42 . To facilitate this switch, the covariance matrix P is recorded from 
the batch estimate 


P k = (. A T k W (23) 

where for the first sequential estimation step the k subscripts are simply the matrix values 
from the batch estimate. For subsequent steps the k subscripts will indicate values from 
the previous step while k + 1 will indicate updated values. A linear Kalman filter is then 
employed to update the estimates of C. As new density measurements are made available 
at one second intervals the estimate of £ is updated using 

&k+ 1 = + 1 + ; + A k + ]PkA% + j) {Y^ + j - A k + j£ k } . (24) 

For the sequential estimator W is just the scalar value of w { calculated using Eq. (19). A is 
only the new row of the A matrix shown in Eq. (16) calculated using the current value of h. 
Y is the current normalized density measurement. To prepare for the next iteration the co- 
variance matrix is updated using 
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P k + 1 ~ P k P k A k + l^k + l +A k + ] P k A k + l^ A k + l P k * ( 25 ) 

This process of updating the polynomials of the density estimate and then updating 
the covariance matrix in preparation for the next step is repeated at one second intervals 
until atmospheric exit. With the density estimate now in place, an estimate of the velocity 
loss to occur in the exit phase due to aerodynamic drag may be computed by integrating 
the drag equation. Begin by writing the drag equation. 


dV ' = J v 2 SC 0 = J v2 r 

dt 2 P r m 2 M d 


(26) 


where M D is the vehicle ballistic coefficient 



m 

c^s- 


(27) 


Replace dt with — in Eq. (26) and substitute the expression given in Eq. (12) for p 
h 

and rearrange terms to obtain 


dv r 0.5 p(h) 2 6 

— =- = — (cj + c 2 x + c,x z + ... +c 7 x°)dh. (28) 

V; M D h 


Since x = h/h £ , if h is a function of h, this expression can be integrated analytical- 
ly between any two altitudes to determine the change in velocity due to aerodynamic drag. 
Again, as was done in the APC algorithm and in Eq. (7), h is chosen to be a constant. If h 
is less than 250,000 ft, p has a discontinuity at = 250, 000 ft; so, the integration must 
be performed in two steps, one from the current altitude to 250,000 ft and a second from 
250,000 to the exit altitude. If we change the variable of integration from h to x, we get 


0.5 p ; 
M D h 




\X S 


-(*-*;) 


r h 


e 

hSl. 


(Cj + CjX + CjX 2 + ... +c 7 x 6 )dx (29) 
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f 1 , . 

[ h el 

\ 

+ 

, -<*-*/) 

) e 

LhS2j 

(Cy + C 2 * + CjX 2 + ... + CyX 6 ) dx 


Kx, 


) 


where x } = h } /h e and the current altitude is expressed x i = h/h g . When h is greater 
than 250,000 ft this integration may be carried out in a single integration step. 


1 0.5 Pj 


M n h 

x D 


\ 


x i 


rM 


( Cj + c 2 x + CjX 2 + ... + c 7 x 6 ) dx 


(30) 


To obtain density from this expression, integration by parts is carried out repeatedly 
to obtain 


jpdx = - 


p.e 


-(X-Xj) 


hS 


r rwi 


-hS~ 

2 

~hS~ 

3 

~hS~ 

5 

^3 

M 

u + 

_ h e . 

u' + 

A. 

u" + ... + 

. h e . 

n 

H 

^3 


(31) 


where 


u = Cj + CjX + c jX 2 + c 4 x 3 + c $x 4 + c 6 x 5 + c 7 x 6 (32) 

and the primes indicate a derivative with respect to x so that 

u' = c 2 + 2c 3 x + 3c 4 x 2 + 4c$x 3 + 5c 6 x 4 + 6c 7 x 5 (33) 

u" = 2c 3 + 6c 4 x + 12c^x 2 + 20c 6 x 3 + 30c 7 x 4 (34) 

This process is continued until all six of the required derivatives are formed using the val- 
ues of c which were most recently estimated, u and all six of its derivatives are calculated 
for the current altitude, and the atmospheric exit altitude. Additionally, u and the six de- 
rivatives must be calculated at 250,000 ft altitude if the current altitude is below 250,000 
ft. These values are inserted into Eq. (31), which is in tum inserted into Eq. (29) or (30) as 
required. The predicted velocity loss due to aerodynamic drag is then found by solving 
Eq. (29) or (30) for the predicted exit relative velocity and then subtracting the current rel- 
ative velocity 
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A^= V„-V r 


(35) 


If the current altitude is below 250,000 ft 


V rx = 


; 0-5pj 

V r + ~M^ 




hSl 


u + 


hSll 2 , 

— — U + ... 4* 
h e\ 


rhSJ l 7 d 6 u 


e J 


dx 6 \ 


(36) 
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or if the current altitude is above 250,000 ft 

( 
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hS2 r hS2 


L h e 


u + 


' e J 


u' + ... + 
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•(37) 


ijr . 


Improved Exit Phase 


This improved exit phase, first published by the Charles Stark Draper Laborato- 
ry 43,44 , is a simplified method of calculating the required altitude rate h for the APC con- 
troller. It is intended to replace the exit phase for the APC algorithm with a simpler 
calculation. To begin the development, the velocity loss due to aerodynamic drag is calcu- 
lated starting with the differential equation for drag: 


dV _ q 
dt M d 


(38) 


Rearrange terms in Eq. (38) to obtain 


dV 


-J-dt. 


M 


D 


(39) 


Replace dt with ^ and expand q. 
dh 
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dV 

V 2 



M D h 


(40) 


With the assumptions of a constant altitude rate and an exponential atmosphere of 
known scale height the above equation may be integrated analytically to obtain the change 
in velocity AV which will occur due to aerodynamic drag. This result is slightly different 
from the original APC exit phase derivation which uses this equation to predict the veloci- 
ty at exit instead of computing the change in velocity which will occur due to drag. The 
preferred form for the drag equation is 


AV = -1/ 


f h des M D 

hSq 



(41) 


To use the hybrid density estimator replace the expression for AV given in Eq. (41) 
with the expression given in Eq. (10). Likewise, to use the polynomial density estimator 
replace the results of Eq. (41) with those of Eq. (35). 

The desired velocity for a vehicle in a purely Keplerian (no aerodynamic forces) or- 
bit at the current radius with the desired altitude rate h to attain the targeted apocenter ra- 
dius may be computed 


des 


^ ^target 

^des 


R 2 

, target , . 

V R ’ 


(42) 


The first term under the radical is the velocity at pericenter for an elliptical orbit with peri- 
center radius R and apocenter at R, arget 


V 

per 


2\lR 


target 


I R ( R + R target) 


(43) 
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A new variable if aaor was introduced 


0.5 


'factor 


' R 2 \ 

<-*££) -i 


per 


( 44 ) 


Therefore V^ es may be written 


des 





(45) 


To avoid the square root in Eq. (45) a small term is added under the radical to complete the 
square 


(r 


des 


* V 


Per, 


+ 2 


' factor V; 2 


v per 


r factor 

V 


r 2 


(46) 


Vj may now be approximated 


V, m V 
des per 


1 + 


9 factor 

“V . 


^5 


+ 0 ( £/ ) = V + ffacto/des + < 47 > 


The corrector step to update altitude rate is a single step Newton iteration. The dif- 
ference between the current inertial velocity minus the velocity loss expected from aero- 
dynamic drag and the desired velocity computed above is called the velocity miss or 

V • . 
miss 

V n,iss = - (V per + r, ac ,y d ' S ) (48) 

dV . 

The negative of V • is then divided by — to produce an update for hj es . 

dhde$ 
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hdes update ^d( 


. (V + / /acwr4 5 )"( V /- AV > 

“des ? 


m d av 2 


— A • 

Two Zr factor r des 


(49) 


Equilibrium Glide to Exit Phase Transition Velocity 

To minimize total AV required to transition from the intermediate orbit to the desired 
orbit it is sufficient to minimize the exit flight path angle provided the vehicle exits in the 
desired orbit plane and the apocenter of the intermediate orbit equals the desired apocenter. 
This approach will maximize the pericenter of the post-aero braking orbit. If the controller 
is able to properly target the apocenter altitude, then minimizing y x will produce a maxi- 
mum exit velocity, a maximum pericenter for the intermediate orbit and a minimum AV. 
Fig. 1 1 shows how selecting a higher transition velocity for the APC controller to switch to 
the exit phase control algorithm will tend to minimize y and the A V required to attain the 
desired final orbit provided the vehicle can properly target the desired apocenter. When the 
transition velocity is increased the predictor/corrector step will calculate a lower h to target 
the desired apocenter. The drawback to minimizing y x by increasing the transition velocity 
and using a shallower flight path for the exit phase is that by doing so the exit phase will be 
flown using a higher percentage of the available lift to follow the desired trajectory. In the 
limit the minimum A V path flies the entire exit trajectory with a bank angle of 1 80° . When 
the transition velocity becomes too great the vehicle can no longer maintain the required 
shallow flight path, even in a relatively smooth atmosphere, and may overshoot the desired 
apocenter altitude, as seen in Fig. 12. Following this type of shallow flight path angle tra- 
jectory severely limits the robustness to density dispersions. If in the initial phases of the 
exit phase, the control system calculates a shallow exit trajectory, one which requires al- 
most full lift down to maintain, any decrease in atmospheric density from that modeled in 
the predictor step will result in less control authority and an inability to fly the shallow tra- 
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Fig. 11 Exit Flight Path Angie and AV Required vs. Transition Velocity 


jectory, less velocity loss than predicted resulting in a faster exit speed than desired and a 
post-aerobraking apocenter higher than desired. An increase in AV results. On the other 
hand, transitioning to the exit phase at a velocity which is too slow guarantees an increase 
in AV by requiring a steep ft to target apocenter which produces large exit flight path an- 
gles. The best trajectory is one which strikes a desirable balance between minimizing A V 
while retaining enough control to be robust under the influence of off-nominal density vari- 
ations. It would seem to be a simple matter to pick a transition velocity which produces the 
desired balance, but the “correct” transition velocity varies with the state of the atmosphere, 
the initial conditions, and the vehicle configuration. 
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Fig. 12 Apocenter and Pericenter Altitudes vs. Transition Velocity 


The two most important parameters used to select an appropriate transition velocity 
are the drag coefficient of the vehicle and the atmosphere yet to be traversed. These two 
parameters and h define the velocity loss that will occur. After considerable testing a de- 
sired altitude rate, h c , of 450 ft/sec was found to yield a good trade off between minimiz- 
ing AV and producing robustness. The simulations depicted in Fig. 1 1 and Fig. 12 require 
a transition velocity of 14,922 ft/sec to produce an exit phase altitude rate of 450 ft/sec. 
As seen in Fig. 11, this transition velocity (and hence this altitude rate) are slightly re- 
moved from the region where exit flight path angle and AV increase dramatically. Yet this 
altitude rate was still steep enough to provide a measure of robustness against density vari- 
ations. Armed with this choice for altitude rate, a better way to calculate transition veloc- 
ity may be formulated. 
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The desired velocity Vj es for a vehicle in a purely Keplerian (no aerodynamic 
forces) orbit at the current radius with the desired altitude rate h to attain the targeted apo- 
center radius was computed in Eq. (47). By adding the velocity loss expected due to aero- 
dynamic drag the current velocity required to target the desired apocenter altitude by 
following a 450 ft/sec path may be computed. The chosen method of density estimation 
may be used to compute the velocity loss due to aerodynamic drag, by inserting the de- 
sired 450 fit/sec altitude rate into the appropriate derivation. Equation (10) is used if the 
hybrid density estimator is the selected method of density estimation, whereas, Eq. (35) is 
used for the polynomial density estimator or Eq. (41) for the simple estimate of a constant 
scale height exponential atmosphere. One additional term is added to allow for the veloc- 
ity loss between initiation of the exit phase and achievement of the desired altitude rate. 
The appropriate velocity to transition from the equilibrium glide phase to the exit phase 
may now be expressed 

^ trig / ~ ^des + ^ (drag) + ( 50 ) 

5r in this equation is the time required from initiation of the exit phase until the de- 
sired altitude rate is attained. The vehicle modeled in this study has a limit of 5°/sec 2 on 
roll acceleration and 20°/sec on roll rate. A value of 20 seconds was selected for 8r be- 
cause with these current limits on roll rate and roll acceleration the vehicle requires thir- 
teen seconds to perform a 180° rest to rest maneuver. After rolling to the lift up 
configuration there is still an additional delay of five to ten seconds before the vehicle’s al- 
titude rate matches the desired value. With 6/ set to 20 seconds the transition velocity cal- 
culation performed extremely well. The methodology for employing a variable transition 
velocity is to compute V (rig/ using eq. (50). When the inertial velocity decreases below 
the calculated W the controller initiates the exit phase. 
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CHAPTER III 

LYAPUNOV CONTROLLERS 

Lee and Grantham present a Lyapunov Steepest Descent controller 7 which is robust 
to at least some atmospheric perturbations. Their controller is for a vehicle which modu- 
lates angle of attack while the MRSR vehicle under study flies at a constant angle of attack 
and varies the bank angle to control the trajectory. A similar controller is developed to 
control the MRSR vehicle. A desired target state is defined for the vehicle at atmospheric 
exit which will minimize the A V required to transition to the desired final orbit. A positive 
definite Lyapunov function is defined such that the vehicle’s state is at the target when the 
Lyapunov function is zero. The control variable is then selected so that the Lyapunov 
function is driven, in a steepest descent fashion, toward the origin. When this method 
failed to be as robust as hoped, a new Lyapunov Tracking controller was developed. 

The Lyapunov Tracking controller permits the introduction of a preferred path lead- 
ing the vehicle to an exit state which gives an acceptable AV to transition to the desired fi- 
nal orbit. In the particular case studied here, the preferred path is recomputed for each 
trajectory based on accelerometer data fed back to the controller early in the flight and a 
“best guess” of the density function for the remainder of the trajectory. Again, a positive 
definite Lyapunov function is defined such that, if the vehicle is on the preferred path, the 
Lyapunov function is zero. The control variable is again selected in a “Lyapunov Opti- 
mal” fashion to drive the Lyapunov function toward the origin as quickly as possible. A 
gain scheduling scheme defines an optimal descent function for each phase of the trajecto- 
ry. Finally, because of high trajectory loads generated by this control scheme and difficul- 
ty in acquiring the desired path, this Lyapunov tracking controller was employed only 
during the exit phase following the modified equilibrium glide 9, 34 phase of the MPC con- 
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troller presented in Chapter 2. The equilibrium glide phase was developed to minimize 
trajectory loads and is very good at doing just that. With the modifications suggested in 
Chapter 2 the equilibrium glide phase control algorithm is very robust to perturbations in 
atmospheric density. The transition velocity from the equilibrium glide phase to this exit 
phase is chosen using the methods presented in Chapter 2 so the trajectory is at the base of 
the preferred path when transition occurs. The Lyapunov Tracking exit phase then follows 
the computed path to exit the atmosphere with exit state very near the minimum A V exit 
state. 


Lyapunov Steepest Descent Controller 


Equations of Motion 

Derivation of the Lyapunov Steepest Descent control algorithm begins with the 
equations of motion for planar flight 

— = ^ = Vsiny (51) 

dt dt 


dV 

dt 


-c D psvl 

2m 



(52) 


dy 

dt 


-C L pSV 2 r 
2 mV 


ll V 

cosO- (-^r--)cosy 
Vr r 


(53) 


Eq. (51) is simply the radial velocity in terms of the inertial velocity and flight path 
angle. Eq. (52) gives the time rate of change in velocity composed of two parts: 1 ) the ve- 
locity loss due to aerodynamic drag and 2) the change in velocity due to gravitational ac- 
celeration, often referred to as the inertial component. Similarly, Eq. (53) is the time rate 
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of change in the flight path angle, also composed of two parts: 1) the change in flight path 
angle due to the component of aerodynamic lift in the vertical plane and 2) the change in 
flight path angle due to gravitational acceleration (the inertial component). The control 
variable O, the bank angle, determines the amount of lift exerted in the vertical plane to 
bend the trajectory and change the flight path angle. 


Nondimensional State Variables 

Dimensionless state variables are introduced: 



X 1 


h/h e 

X = 

x 2 

= 

v/{Jp/R) 


X l 


7 


(54) 


along with a dimensionless time variable x 

x = (t/h e ) Jvl/R. 

The equations of motion may now be written: 

x j = A^sinx^ 


(55) 


(56) 



r sinx ? 

(c- 1 +Xj ) 2 


(57) 


Aaxi 


x 3 = 


-cosO + 


COSJtj 

c- 1 + X, 


X 2 - 


( C-1 +Xj)x 2 ] 


(58) 


where a = p/p 0 = exp [(- (h - h Q )) /hS] , A - (p 0 Sh e C L ) / (2m) , 

B = (P 0 Sh e C D )/(2m) and c = R/h £ . It has been found that a good approximation is 
to assume that the relative and inertial velocity differ by a constant, so that 
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and similarly 


V r = V-5V 


(59) 


X 2 r ~ x 2~& x 2 - (60) 

For this controller it is more convenient to replace cos<t> in Eq. (58) with the control 

variable u where u = cosO and is thus bounded between ±1. Eq. (58) is therefore re- 
placed with 


A OX 2 u 

X, = — + 


COSJTj 

c-J+x 


1 L 


x 2 ~ 


(c - 1 + Xj) X 2 


( 61 ) 


Target State 


The minimum AV aerobraking maneuver is one which exits the atmosphere on a tra- 
jectory with the correct apocenter altitude and a maximum vacuum pericenter altitude. 
This goal is attained by exiting the atmosphere with the minimum possible flight path an- 
gle and the correct velocity to attain the desired apocenter. The goal, therefore, is to guide 
the vehicle along an aerobraking trajectory which reaches the atmospheric interface alti- 
tude with the correct velocity to attain the desired apocenter altitude while maintaining a 
minimum positive flight path angle at exit. The flight path angle must remain positive for 
the vehicle to exit the atmosphere. This design objective is established by setting the tar- 
geted flight path angle at atmospheric exit to zero and establishing a target exit velocity. 
The target state may be presented in non-dimensional form as 7 


x = 



(62) 


The target exit velocity, and hence x 2 may be derived assuming a Keplerian orbit 
from atmospheric exit to apocenter. This desired exit velocity is a function of the exit 
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flight path angle, and several constants for the problem including the atmospheric interface 
radius (/?) and target apocenter radius (r fl ) . The desired exit velocity is 


*2 ~ 



r r 1 

f 

(< 2 « 

['-*] 

n 




(63) 


Descent Function 


A function is a descent function if, and only if, it is a positive definite differentiable 
function. That is: 4 ^ 


W (*) >0 for all x * X 

(64) 

W{&) = 0 

(65) 

W(x) 

— 5 — * 0 for all x * x 

OX 

(66) 


Any candidate Lyapunov function may be chosen as the descent function W [x (/) ] . 

n 

However, the most logical choice, and the one recommended by Lee and Grantham , is a 
weighted quadratic measure of distance to the target. This function is expressed 


W(x) = 

xj-l x 2 -x 2 * 5 ] 

Pll °P]2 

0 10 

Xj-l 

x 2 -x 2 



P 12 °P 2 2_ 

_ X 3 _ 


where the constant weighting terms p ^ are chosen to define a preferred direction toward 
the target in the - Xj state space. The preferred direction for the states is presented in 
Fig. 13. An ellipsoid is chosen, oriented so that, while the vehicle is deep in the atmo- 
sphere, the preferred direction (opposite the descent function gradient) in the Xj - x 3 state 
space gives positive lift to climb out of the atmosphere, but as the vehicle approaches at- 
mospheric exit the preferred direction uses negative lift to minimize the exit flight path an- 
gle. The weights must be scaled so the velocity reaches the target velocity as the vehicle 



Fig. 13 Preferred Xj - x 3 Direction of Motion 
(Adapted from Reference 7) 


reaches the atmospheric interface altitude. For the elliptical descent function shown in 
Fig. 13, p xx may be calculated as follows 7 . 

2 2 2 2 

Pjj = a sin ty + b cos (68) 

p 12 = sin <()cos (a -b ) (69) 

Pll = b 2 sin 2 § + a 2 cos 2 ty (70) 


The angle between the gradient of the descent function and the state space velocity 

vector f(x,u) is expressed 

ri/ x _ (dW(x))/(dx) f(x, u) _ a 

H{x,u) it (aw(jf))/(9oc) ii \\fix, u) ii cosP (71) 

where P is the angle between the gradient of the descent function and f(x, u ) . 
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Lyapunov Steepest Descent Optimal Control 

♦ 7 

For the control to be Lyapunov steepest descent optimal, u ( x ) must make 

H [x, u (x)]<H [x, u (x) ] (72) 

for all u e U where U is the allowable set of controls bounded by ±1. Furthermore, for 
the control to be Lyapunov steepest descent optimal, f[x, u* (x) ] * 0. If it were possible 
to make H [x, u (x) ] <0 everywhere, then global stability with respect to the target 
could be guaranteed. Even if H [x, u* (x) ] <0, asymptotic stability with respect to the 
target could be guaranteed. Unfortunately, with u bounded between ±1 neither of these 
outcomes is always possible. Even so, u * (x) tries to move the system state variables as 
nearly opposite the gradient of the descent function as possible, given the dynamics of the 
system and the limits on the control. 

To determine u (x) set =0 and solve for u. If this value of u lies between ±1, 

du 

then u* (x) is either this value or ±1, whichever minimizes H. If the value of u which 

solves ^ = 0 is not between ±1, then u (x) is selected from±l to minimize H 47 . 

du 

Performance Results 

The Lyapunov Steepest Descent feedback control algorithm will guide the vehicle 
to very near the minimum AV exit state provided the weights, and hence, W (x) is 
properly selected. Unfortunately, these p ^ weights must be readjusted to attain accept- 
able performance for each perturbed entry condition, vehicle lift and drag perturbation, or 
atmospheric density perturbation. No acceptable method, other than a manual search, was 
found to determine the appropriate weighting for each perturbed run. Clearly, this lack of 
asymptotic stability is not compatible with the objectives of this research. 
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An appropriate descent function is found for this controller to perform acceptably for 
the nominal case of a vehicle targeting a 270 nm circular orbit after entering the Martian 
standard atmosphere with 6.0 km/sec velocity relative to the planet, -12° entry flight path 
angle, and a lift to drag ratio of 1 . The same vehicle and the same entry conditions are also 
simulated, assuming both a high and a low density Martian atmosphere. The ellipse which 
determines the descent function is chosen to have a semimajor axis of 1 .65, a semiminor 
axis of 0.41, and a rotation angle <]) of 4.2°, measured as shown in Fig. 13. A few pertur- 
bations from the nominal case are then simulated and the somewhat disastrous results are 
presented in Table 1 along with the optimal results for the same perturbations generated 
using the method of Appendix A. 

This controller is simply not very robust, given density, navigation or vehicle pertur- 
bations expected for the Martian aerobraking problem. The controller may be fine tuned 
for one rate of energy depletion, but if anything alters the rate of energy loss the controller 
must be readjusted, by altering the relative weights between the states, to bring the veloci- 
ty to the targeted velocity exactly as the vehicle passes through the atmospheric interface 
altitude. A steeper entry flight path angle thrusts the vehicle deeper into the atmosphere, 
thereby increasing the rate of energy loss. Likewise, an atmosphere which is more dense 
than expected, or a drag coefficient higher than expected, causes the vehicle to lose energy 
at a higher rate than planned. The resulting exit conditions are too slow and at too low an 
apocenter altitude. In the worst trajecories the vehicle fails to exit the atmosphere at all. 
Similarly, a shallower entry flight path angle, less dense atmosphere or lower drag coeffi- 
cient results in less velocity loss than intended and apocenter altitudes higher than desired. 
To reduce the sensitivities to perturbations which change the rate of energy loss the Ly- 
apunov controller is therefore reformulated as a tracking controller designed to follow a 
chosen path to atmospheric exit. 



Table 1 Lyapunov Steepest Descent Controller Performance Results 
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Lyapunov Tracking Controller 

To gain acceptable robustness for this controller, the methodology is changed from a 
steepest descent controller which targets the optimal terminal state to a steepest descent 
controller which targets a preferred path. That path then is selected to lead the vehicle to a 
desirable exit state with enough robustness to prevent minor density upsets from being cat- 
astrophic. 

The Preferred Path 

Derivation of this controller begins with definition of the preferred path. As with the 
predictor corrector algorithms a constant altitude rate path leading to the desired atmo- 
spheric exit state is selected. The difference between this Lyapunov Tracking Controller 
(LTC) and the predictor correctors is in how the controller computes the constant altitude 
rate path. The predictor corrector algorithms use various methods to select a constant alti- 
tude rate which will give the desired apocenter altitude and then use altitude rate error to 
select the appropriate bank angle. The LTC, on the other hand, assumes that it is desirable 
to always fly the same altitude rate to atmospheric exit and arrive there with the appropri- 
ate velocity to achieve the proper apocenter altitude. The LTC then chooses the in-plane 
portion of lift to approach the path in a steepest descent fashion. The altitude rate is select- 
ed to produce the desired trade-off between robustness to density perturbations while still 
minimizing the total AV required. 

A constant altitude rate of 450 ft/sec is again selected (as on page 37) to define the 
desired path leading to atmospheric exit with the appropriate velocity to target the desired 
apocenter altitude. This altitude rate produces trajectories which require within 20 to 30 
ft/sec of the minimum A V values for the various expected perturbations without short pe- 
riod density upsets, yet is still robust to density variations of ±50% over small altitude in- 
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tervals. The equations used to derive the improved exit phase for the Mars Predictor 
Corrector are employed again here to define this path. The velocity required at a given al- 
titude flying a specified altitude rate assuming a Keplerian orbit (no aerodynamic drag ef- 
fects) was given in Eq. (47) but is repeated here for completeness. 


des 




= V + r 


_1_ r> *• 

per factor des 


(73) 


The velocity loss expected due to aerodynamic drag is added to this velocity to deter- 
mine the current velocity for the desired path. Note that this desired velocity is a function 
of the dynamics of the Martian orbit, the current altitude, the selected altitude rate (450 ft/ 
sec) and the expected velocity loss (which is a function of the expected atmospheric densi- 
ty function and the vehicle coefficient of drag). The velocity loss expected due to aerody- 
namic drag is calculated assuming the 450 ft/sec altitude rate path will be flown using Eq. 
(10) of the hybrid density estimator, or with Eq. (35) of the polynomial density estimator, 
or with Eq. (41) using the simplification of a constant scale height exponential atmo- 
sphere. The desired current velocity defining the preferred path is 


V = 




(74) 


This velocity may be converted to non-dimensional form 

*2 = V/{.J)JR) (75) 


The desired flight path angle % 3 is computed 

x 3 = asin(h des /V). (76) 


Together, x 2 and x 3 define the preferred path which leads the vehicle along a robust 
corridor to a desirable exit state. Now, a Lyapunov function must be formulated and a 
control found which will drive the vehicle onto and then along the chosen path. 
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The Lyapunov Descent Function 


The selected positive definite Lyapunov function is 


W(x) = C 

c 2“*2 x 3-*3\ 

P 11 P 12 

X 2 -X 2 

L 


P 12 P22_ 

x 3~ X 3_ 


( 77 ) 


This function is analogous to distance from the target path and is zero whenever the vehi- 
cle is on the target path and positive otherwise. Again, the p ^ values are chosen to form 
an ellipsoid, the negative gradient of which defines the preferred approach to the target 
path. This ellipsoid is shown in Fig. 14. The p xx values are computed from the semima- 



Fig. 14 Lyapunov Tracking Controller x 2 - x 3 Descent Function 
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jor axis a, the semiminor axis b y and the angle of rotation <)> of the ellipse defining the pre- 
ferred gradient 7 onto the chosen path as in Eqs. (68) through (70). They are repeated here 
for completeness 

p }] = a 2 sin 2 4> + b 2 cos 2 4> (78) 

Pj 2 - sin§cos§(a 2 — b 2 ) (79) 

p 22 - b 2 sin 2 <J> + a 2 cos 2 <J> (80) 


Again, the angle between the gradient of the descent function and f{x, u) is ex- 
pressed 


H(x,u) 


(3W(jc) ) / (3 jc) f(x,u) 
(dW(x))/(dx) || || f(x, u) || 


(81) 


where P is the angle between the gradient of the descent function and the state space ve- 
locity vector /(x, u ) . 


Selection of the Control 

As was done earlier in the report for the Lyapunov Steepest Descent Controller, a 
control it (x) is sought which will move the system state variables as nearly opposite the 
gradient of the descent function as possible, given the dynamics of the system and the lim- 
its on the control. 

As before, to determine it (x) set ^ = 0 and solve for u. Note however, that H in 

this discussion is not the same function as H in the LSDC discussion. If the value of u lies 

between ±1 then u (x) is either this value or ±1, whichever minimizes H. If the value of 

u which solves ^ = 0 is not between ±1, then it (x) is selected from ±1 to minimize H. 
du 
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Performance Results 

An acceptable descent function for this algorithm has not been found. Fig. 15 
shows the x 2 and x 2 time histories while Fig. 16 shows the x 3 and x 3 time histories for a 
trajectory guided by this LTC. The elliptical descent function used in these simulations 
has a semimajor axis of 40, a semiminor axis of 1, and a rotation angle <|> of 19°, measured 
as shown in Fig. 14. This descent function produces an apocenter altitude following the 
aerobraking maneuver of 356 nautical miles which is as close to the 270 nm target as 
could be obtained while keeping the flight path angle near the optimal. But as Fig. 15 
shows, as the vehicle neared the target path, the LTC did not make the final correction nec- 
essary to converge on the desired velocity. A second descent function is formed which 
guides the vehicle closer to the target apocenter altitude. This descent function uses an el- 
liptical function with the same semi -major axis of 40, a semi-minor axis of 1, but the rota- 
tion angle is changed to 55°. As Fig. 17 and Fig. 18 plainly show, the desired apocenter 
altitude is not achieved by following the desired path to exit but rather by reducing the ve- 
locity ( x 2 ) more than desired and then climbing with a steeper flight path angle ( x 3 ) than 
preferred. Almost by accident, the desired apocenter altitude is attained in this one case. 
No fixed configuration for the elliptical descent function was found which consistently al- 
lows this algorithm to acquire the target path and follow it to an acceptable exit state. 

Intuitively, it is easy to see that, if the velocity is slower than V and the flight path 
angle is less than % 3 , the logical choice is to use positive lift to get closer to the path. Like- 
wise, if V is greater than 9 and the flight path angle is greater than % 3 a lift down orienta- 
tion is required to approach the path. The ambiguous areas are in the other two quadrants, 
where either the velocity is too fast but the flight path angle is too shallow, or where the 
velocity is too slow but the flight path angle is greater than desired. It is desirable to define 
a line passing through the target state at each instant in time separating the regions. The 
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slightly different altitude rate. Though the 450 ft/sec altitude rate that was chosen gives 
excellent results, it should still be acceptable to fly a 425 or a 475 ft/sec altitude rate path 
which would lead the vehicle to the desired apocenter orbit. 


To define this switching line, the angle of rotation <(> (Fig. 14) of the descent func- 
tion is varied during the trajectory such that 


<)) = atan 


r d & 3 

dh_ 

dx 2 

L dh J 


(82) 


In effect this expression defines a switching line formed by linearizing about the current 
target state and varying altitude rate. Though 


dW(x, t) 

dx 


... ,, dW 

(f(x, u)) , 


(83) 


dWd<b dW 

because of the missing ^ component. But, since ^ is small compared to the ele- 

dW ^ V 

ments of and since <}> varies slowly, increasing monotonically from about 15° to about 
ox 

75° during the exit phase, this component is assumed to be insignificant. The commanded 
bank angle is still determined as before by selecting the value of u which minimizes H(x,u) 
with H(x,u) defined as Eq. (81). 


Though this method of varying the weighting matrix gives improved performance, 
the algorithm still has problems acquiring the target path. The vehicle still uses lift down 
too early and plunges deeply into the atmosphere, creating extremely high vehicle acceler- 
ations and heat rates in the process. To cure this problem the Lyapunov Tracking Algo- 
rithm (LTA) developed here is incorporated as an exit phase following the equilibrium 
glide phase of the MFC algorithm outlined in Chapter 2. 



56 


Lyapunov Tracking Controller Exit Phase 

The equilibrium glide phase was developed to guide the vehicle into the atmosphere and 
hold it in equilibrium until the velocity has been appropriately reduced. It must perform this 
task while keeping the maximum trajectory loads and peak heat rates within bounds. It per- 
forms this task very well. On the other hand the LTC just described does a good job of holding 
the desired path to exit if somehow it is started near that path. If we break up the trajectory so 
that the LTC has responsibility for control when these conditions exist, its strong points are 
likely to improve overall performance. With this idea in mind the marriage of the LTA as an 
exit phase with an equilibrium glide phase controlled with one the controllers developed earli- 
er is a natural implementation. Such a concept is implemented next. The method of comput- 
ing transition velocity from the equilibrium glide phase to the exit phase presented in Chapter 
2 placed the vehicle very near the desired path. This combination of equilibrium glide phase 
and LTA proved to be the best controller examined. 

The two density estimation techniques presented in Chapter 2 are also tested. The com- 
plete control algorithm with the hybrid density estimator included will be referred to as the 
Hybrid Lyapunov Tracking Controller (HLTC) while this controller with the polynomial den- 
sity estimator will be called the Lyapunov Tracking Controller (LTC). 
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CHAPTER IV 

CONTROLLER SENSITIVITY ANALYSIS 

The performance of the MPC and MHPC control algorithms developed in Chapter 2 
and of the LTC and LHTC control algorithms developed Chapter 3 was determined along 
with the performance of the APC and Energy Controllers which are derived and discussed 
in detail in the literature 9 ' 101213,48 . All the algorithms are tested using a six degree of free- 
dom computer simulation based on the Program to Optimize Simulated Trajectories^ 6 
(POST), which utilizes a fourth order Runge-Kutta numerical integration scheme to con- 
tinuously integrate both the force and the moment equations of the vehicle. The control 
algorithms are tested to determine the effect of large scale density variations such as those 
caused by the seasonal sublimation and condensation of the Martian atmosphere or by a 
global dust storm. They are also run to determine the effect on the performance of short 
period atmospheric variations by injecting square wave density pulses, similar to those 
used by Fitzgerald 1 ®’ **, of various magnitudes and durations into the density function at 
various altitudes. Entry flight path angles are varied within the current predicted error 
band 46 . Perturbations in the vehicle lift and drag characteristics are also simulated. Final- 
ly, combinations of these perturbations in the atmospheric density function, entry flight 
path angle and vehicle lift and drag characteristics are examined and the performance of 
each controller is analyzed. 

Following a brief description of the vehicle and trajectory simulation program used, 
the data from this test program are presented graphically utilizing three dimensional mesh 
plots. The primary thrust of this simulation effort is to select the best controller(s) from 
those studied. A fuller performance evaluation of the selected controllers, aimed at deter- 
mining their robustness limits, is presented in Chapter 5. 
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Vehicle and Trajectory Simulation Inputs 

The vehicle used in the study is a biconic aeroshell design with a fifteen foot base di- 
ameter and a weight of 1 1,023 lbs. The base surface area of 176.1 ft 2 is used as the refer- 
ence surface area. The vehicle is designed with a five foot center of gravity offset resulting 
in a trim angle of attack of 27° which is maintained throughout the maneuver via a simple 
proportional feedback controller. Control is attained through bank maneuvers which re- 
orient the direction of the lift vector. These bank maneuvers are commanded as body axis 
rolls with coordinating body axis yaw maneuvers. The nominal lift coefficient is 0.68892 
while the drag coefficient is 0.69819, producing a nominal L/D of 0.99. 

The Mars Global Reference Atmosphere Model 35 (MARS-GRAM) is used to pro- 
duce realistic atmospheres for the study. Three different atmospheres representing a nom- 
inal, a low density and a high density Martian atmosphere are considered (Fig. 19). The 
nominal atmosphere is the COSPARV Model Atmosphere For Mars 18 , while the high den- 
sity and low density atmospheres are derived using MARS-GRAM. The low density at- 
mosphere is a MARS-GRAM simulation of the lowest density Martian atmosphere 
predicted for April 10, 1999 assuming no dust storms and a 10.7 cm solar flux of 50 (nom- 
inal value = 150). The high density atmosphere represents the highest density atmosphere 
predicted on December 27, 1997, again with no dust storms but this time with a 10.7 cm 
solar flux of 300. Although MARS-GRAM was originally incorporated as a subroutine to 
POST which could be called to generate atmospheric data on line, MARS-GRAM is not 
utilized in this manner because of the added computational time. MARS-GRAM gener- 
ates atmospheric data which is stored in tabular form. These tables of atmospheric data 
are then included in the POST input namelist. 

In addition to the large scale density variations introduced by using the low, nominal 
or high density atmosphere models described above, short period variations in the atmo- 
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Fig. 19 Mars Nominal, Low and High Density Atmospheres 


spheric density function are investigated by introducing square wave density pulses, simi- 
lar to those used by Fitzgerald 10, n . These pulses perturb the local atmosphere within a 
10,000 or 20,000 ft altitude band by multiplying the expected density by a constant magni- 
tude density multiplier. The magnitudes of the density multipliers used include 0.5, 0.75, 
1.25, 1.5, 1.75 and 2.0. The lower edge of the density pulses are varied in 10,000 ft steps 
from 100,000 to 290,000 ft 

The atmospheric interface altitude is taken to be 125 km (410,105 ft) and the initial 
conditions are defined at this altitude. The entry velocity is 6 km/sec (19,685 ft/sec) and 
the nominal entry flight path angle is -12°. The targeted orbit is a 270 nm circular orbit. 
In addition to the atmospheric perturbations mentioned above, perturbations are intro- 
duced in the vehicle lift and drag coefficients representing variations of ±33% from the 
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nominal L/D ratio of 0.99. The +33% L/D perturbation is introduced by multiplying the 
nominal lift coefficient by 1.14, while the drag coefficient is multiplied by 0.86. The 
-33% L/D perturbation is introduced by multiplying the nominal lift coefficient by 0.8, 
while the drag coefficient is multiplied by 1 .2. This method of varying UD also perturbed 
the ballistic coefficient of the vehicle. Navigation errors in the form of variations in the 
entry flight-path-angle of ±0.25° and ±0.5° from the nominal -12° are considered. The 
performance for each pertuibed run is presented as total AV required to achieve the desired 
final orbit. AV is a measure of the controllers overall success in meeting the desired exit 
conditions. The AVis calculated assuming one bum at atmospheric exit oriented along the 
velocity vector to correct any apocenter error, a second at apocenter to raise pericenter, and 
a final bum to correct any orbit plane error. 


Analytic Predictor Corrector Performance Results 

The original APC controller 9 - 10 - ’2. 13.48 did not fare very we „ when chaJ]enged whh 
the postulated perturbations used in this study. Figures 20, 21, 22, and 23 present the per- 
formance of this controller when faced with these perturbations. These charts summarize 
the performance of a nominal vehicle which enters the atmosphere with an v = -12° 

*e 

and then encounters a square wave density pulse. The AV required to circularize is plotted 
on the vertical axis. Figure 20 presents the results of encountering square wave density 
pulses in a nominal Martian atmosphere, while Fig. 21 shows similar data for a low densi- 
ty atmosphere and Fig. 22 then illustrates the results for a high density atmosphere. In the 
first diagram of each figure the density pulse perturbs a 10,000 ft altitude band while in the 
second diagram the pulse affects a 20,000 ft band. Magnitudes for these pulses range from 
-50% to +100% in 25% increments. The location of the lower edge of the pulse was 
moved from 100,000 ft to 290,000 ft in 10,000 ft increments. The magnitude of the pulse 
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and the altitude of the lower edge are shown on the diagrams. Figure 23 shows the AV re- 
quired to circularize as LTD and *1 e are varied. Figure 23a shows the results of these varia- 
tions in a nominal atmosphere, while Fig. 23b presents the same perturbations in a low 
density atmosphere and Fig. 23c in a high density atmosphere. 

Figure 23 indicates that the APC controller exhibits a considerable sensitivity to off 
nominal vehicle aerodynamics and to navigation errors. This controller also shows a 
marked decrease in performance in the high density atmosphere with no density steps 
when compared to its performance in the low density and nominal atmospheres. The AV 
required to circularize following the aerobraking maneuver in a high density atmosphere 
even with no density step (0% magnitude density step) is over 400 ft/sec while the optimal 
results presented in Table 1 show that it should require less AV to circularize after aero- 
braking in a high density atmosphere (optimally about 316 ft/sec) than in a nominal or low 
density atmosphere. Additionally, the variations in AV shown in Fig. 23c are considerably 
worse than those in Fig. 23a or b. Part of this sensitivity comes from using a specified 
transition velocity to switch to the exit phase, ignoring the actual energy loss to occur dur- 
ing the exit phase. The other reason for this sensitivity is a rather simplistic density model. 
Also, the APC controller is less sensitive to density steps in the high density atmosphere 
than in the low or nominal atmosphere. This sensitivity can again be explained by the 
choice of a specified transition velocity for the switch from entry to exit phase. 

The transition velocity selected for this controller was 14,922 ft/sec. This transition 
velocity is appropriate for a nominal vehicle which enters the nominal atmosphere with a 
nominal flight path angle of -12°. Also, if the initial flight path angle is steeper than -12° 
or the atmosphere is more dense than expected, the vehicle will plunge into the atmo- 
sphere deeper than expected, and consequently, will have more atmosphere to traverse 
during the exit phase and will lose more energy to aerodynamic drag. Similarly, if the ve- 
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hide’s drag coefficient is higher than estimated, the vehicle will lose more energy during 
exit than planned. Forcing the vehicle to decelerate to a predetermined velocity before ini- 
tiating the exit phase requires the exit phase to be flown at a larger altitude rate than de- 
sired to reach the nominal apocenter altitude. This approach results in higher AV values 
and tends to amplify the effect of off nominal entry conditions or drag coefficient. The 
steeper exit path flown by this controller in the high density atmospheres is, however, more 
robust to density variations, as is suggested in Fig. 22. This result is due to the fact that a 
trajectory which flies a steeper exit phase has reduced more of the velocity deep in the at- 
mosphere and is not requiring as much velocity loss during the exit phase. Density varia- 
tions which perturb the amount of velocity loss during the exit phase have less effect when 
more of the velocity is reduced deep in the atmosphere. Later in the report the value of 
making this transition velocity an adaptive parameter will be shown. 

The second concern with the APC controller is the density estimation technique. 
The density estimator built into this algorithm assumes the density function is a fixed scale 
height exponential function. The density derived onboard from accelerometer measure- 
ments is filtered using a low pass filter to remove high frequency noise. The result is then 
used to bias the exponential function used to estimate density. This technique works well 
as long as the density function does not vary much from a smooth exponential function 
and, more critically, the scale height of the atmosphere is fairly constant and does not vary 
appreciably from the assumed scale height. Unfortunately, the scale height of the Martian 
atmosphere does vary considerably 

Fig. 19 shows the range in the density function predicted by MARS-GRAM. This 
figure presents altitude versus log density. The scale height may be determined by taking 
the negative of the slope of the density function from this graph. The scale height does not 
vary considerably below 250,000 ft, but above 250,000 ft there is considerable variation. 
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This variation does not affect the trajectories flown in the low density atmosphere very 
much because, in a low density atmosphere above 250,000 ft, aerodynamic forces have 
negligible effect on the vehicle. But, for those trajectories flown in the high density atmo- 
sphere, failure to properly model the atmosphere has considerable effect. The estimated 
density falls well short of the actual values above 250,000 ft. Consequently, there is more 
aerodynamic drag than predicted. Typically, this density modeling error leads to a final 
apocenter altitude ten to twenty nautical miles lower than targeted. This problem, com- 
pounded by the steep altitude rate in the exit phase (brought on by the constant transition 
velocity) reduces the controller’s ability to correct for density upsets which occur after the 
exit phase is initiated. A noticeable sensitivity to square wave density pulses (Fig. 20, 
Fig. 21 and Fig. 22) results. 

Overall, the APC is still a good controller. It guides the vehicle through the aero- 
braking maneuver with minimal heat, acceleration, and dynamic pressure loads. It exits 
with a state near the optimal one when the density function encountered is near the nomi- 
nal value, when navigation is good enough to allow precise control over the entry state, 
and when the hypersonic lift and drag characteristics of the vehicle are close to the design 
values. As part of the modification of this controller to meet the Martian requirements the 
value of K- was changed to 4.5 as recommended in Chapter 2. This change kept the vehi- 
cle from exiting the atmosphere before slowing enough to transition to the exit phase 
(skipping out) for all of the test cases examined. However, the APC just is not quite ro- 
bust enough to adequately handle expected perturbations in the Martian atmosphere, vehi- 
cle entry conditions, or vehicle lift and drag variations. 

Energy Controller Performance Results 

Figures 24, 25, 26, and 27 show the performance of the Energy Controller. Again, 
the results are presented in the same format as before with Fig. 24 showing results for den- 
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sity steps in a nominal atmosphere while Fig. 25 shows the results for a low density atmo- 
sphere and Fig. 26 shows the performance results in a high density Martian atmosphere. 
Fig. 27 shows the effect of varying lift-to-drag ratio and entry flight path angle with Fig. 
27a being in a nominal atmosphere, Fig. 27b showing these results in a low density atmo- 
sphere, and Fig. 27c showing the results in a high density atmosphere. 

The Energy Controller is substantially more robust than the APC controller with re- 
spect to vehicle lift and drag variations and to navigational errors. Figure 27 shows practi- 
cally no variation in AV required to circularize. Furthermore, these results all fall below 
400 ft/sec to circularize. This insensitivity may be attributed to the fact that the Energy 
Controller does not assume a density function, though an exponential function is expected; 
instead, it relies on the current energy rate and energy error to determine which path 
should be pursued. Variations in the vehicle’s drag coefficient simply change the energy 
rate and the controller compensates . Likewise, variations in the overall state of the atmo- 
sphere (low, nominal, or high density atmosphere), or variations in the entry flight path an- 
gle which force the vehicle deeper or shallower into the atmosphere are seen by the 
controller as changes in the energy rate. Since the controller seeks to make energy rate ap- 
proach zero as energy error approaches zero, variations of this type are handled well. 

This method works well as long as the density function is a smooth exponential but, 
as the 10,000 and 20,000 ft density pulse diagrams illustrate, the Energy Controller shows 
definite sensitivity to density functions which are not smooth. The large magnitude 
20,000 ft duration density steps are more than this controller can tolerate. Figure 25b 
shows that in a low density atmosphere, 20,000 ft duration density steps of +75 and 
+100% magnitude with lower edges between 100,000 ft and 120,000 ft, cause a cata- 
strophic failure requiring more than 1000 ft/sec of propulsive maneuvering to circularize 
the desired orbit. These failures occur because the vehicle enters the high density region 
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near the bottom of the trajectory. When the onboard accelerometers sense the rapid decel- 
eration caused by the high density pocket and feed this information to the controller, the 
control system responds by applying lift up to decrease the energy rate. As the vehicle ex- 
its the high density pocket the control system responds by commanding a lift down config- 
uration. But the vehicle’s response time is too slow, requiring thirteen seconds to perform 
a 180° rest to rest roll maneuver. By the time the maneuver is complete the vehicle has 
moved higher in the atmosphere and no longer is able to control the trajectory using aero- 
dynamic forces. The vehicle exits the atmosphere without properly depleting the kinetic 
energy. This behavior could also be called a skipout. The same phenomena is observed 
for the high magnitude density pulses in the nominal and high density atmospheres, 
though the effect is less disastrous. 

The locations of the density pulses which cause the problems are higher in the nomi- 
nal atmosphere than in the low density atmosphere, and even higher still in the high densi- 
ty atmosphere than in the nominal atmosphere. These higher critical locations occur 
because the vehicle’s initial configuration is lift up. Higher density atmospheres exert 
more aerodynamic force at higher altitudes, tending to decrease the vehicle’s negative alti- 
tude rate earlier. They also increase the altitude at which the vehicle bottoms out. A den- 
sity pulse which perturbs the trajectory near its minimum altitude must be located higher 
in a high density atmosphere than in a low density atmosphere. 

One additional drawback to the Energy Controller is higher trajectory loads than the 
algorithms which use an equilibrium glide phase. The equilibrium glide phase holds the 
lift up configuration until the trajectory bottoms out in almost all cases. The Energy Con- 
troller will roll the vehicle from lift up before the vehicle bottoms out, allowing the vehicle 
to sink to a lower minimum altitude, producing higher peak aerodynamic heating loads, 
higher maximum dynamic pressures, and higher maximum acceleration loads. These 
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higher trajectory loads may reduce somewhat the advantages of aerobraking, especially if 
they require the vehicle to be more rigid to withstand the acceleration forces, or if they re- 
quire additional heat shields or ablative materials. Though this study is concerned prima- 
rily with control system robustness, the effect of a control system on trajectory loads must 
also be considered. 

The Energy Controller has shortcomings, especially with respect to short period den- 
sity variations and trajectory loads, which make it unsatisfactory for controlling an aero- 
braking vehicle operating in the Martian atmosphere. 

Mars Hybrid Predictor Corrector Performance Results 

The Mars Hybrid Predictor Corrector (MHPC) was one of the two best performing 
algorithms tested for this series of perturbations. As discussed in Chapter 2, this algorithm 
employs a variable transition velocity for the switch from the equilibrium glide phase to 
the predictor corrector exit phase. Equally important is the density estimation technique 
which measures and records density at discrete altitude locations during the entry into the 
atmosphere. Density during the exit from the atmosphere is measured and compared 
against that predicted using the stored entry data. The result is filtered to remove high fre- 
quency noise and used to bias the density estimate developed during entry. This biased 
density estimate is then used to predict velocity loss for the remainder of the trajectory. As 
may be suspected, this method is extremely effective whenever the density profiles for the 
inbound and outbound legs of the trajectory are the same. 

The performance of this controller is presented in Figures 28, 29, 30, and 31 . Again, 
the first three figures summarize the performance when the density function is perturbed 
with 10,000 and 20,000 ft duration square wave density steps. Also, Fig. 28 presents these 
results when perturbing waves are injected into a nominal atmosphere. Figure 29 shows 
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the results in a low density atmosphere and Fig. 30 in a high density atmosphere. In all of 
these cases the controller is able to guide the vehicle to the targeted exit state almost per- 
fectly. With this controller the required AV following a maneuver in a high density atmo- 
sphere is slightly less than that in a nominal atmosphere, which is in turn slightly less than 
that in a low density atmosphere. These results are in agreement with those found using 
the Conjugate Gradient optimization technique outlined in Appendix A. 

The performance when the vehicle lift and drag characteristics are varied and when 
the entry flight path angle is varied (Fig. 31) are equally promising. The reader is cau- 
tioned that the results presented here were all generated with density functions which are 
simply functions of altitude. The density function for the outbound leg of the trajectory is 
identical to the density on the inbound leg. The density estimator in this control algorithm 
gives excellent results when the outbound density function matches that measured while 
inbound, and the control algorithm is able to guide the vehicle to near perfect exit state 
whenever it is supplied with a good density function estimate. Later, in Chapter 5 the per- 
formance will be evaluated whenever the inbound and outbound density functions differ. 

Mars Predictor Corrector Performance Results 

The Mars Predictor Corrector (MPC) Control Algorithm differs from the MHPC of 
the previous section only in the density estimation technique employed. The MPC mea- 
sures and stores density every second during the descent into the atmosphere. These den- 
sity measurements are then normalized using a two stage exponential function. The 
resulting normalized data are fit with a sixth order polynomial in altitude. This polynomi- 
al is continually updated throughout the trajectory after each density measurement is tak- 
en. Again, the resulting density estimate is used to compute the velocity loss yet to occur 


due to aerodynamic drag. 
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The data generated during the simulations of this phase of the evaluation of the MPC 
were not as good as those from the MHPC. This density estimation technique might be 
expected to perform slightly worse, certainly no better, than the hybrid density estimation 
technique whenever the inbound and outbound density functions are the same. The 
strength of this density estimation technique is expected to be for those cases to be studied 
later when the inbound and outbound density functions are different (Again, such simula- 
tions are analyzed in Chapter 5). 

The performance of this algorithm is presented in Figs. 32, 33, 34, and 35. The first 
three of these figures present the results when square wave density pulses are imbedded in 
the Martian atmosphere, while the last figure shows the results of varying lift-to- drag ra- 
tios and the entry flight path angle. 

The performance of this algorithm as depicted in the following four figures, though 
not as good as that of the MHPC, is still quite acceptable. The worst performance noted 
here, caused by a 20,000 ft duration, +75% density pulse in the high density atmosphere 
located between 180,000 and 200,000 ft, required 457 ft/sec to attain the desired orbit. 
This algorithm, when faced with variations in enti7 flight path angle and L/D, produces 
practically flat performance maps that are very near the idealized near-optimal values cal- 
culated using the method of Appendix A. 
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Lyapunov Hybrid Tracking Controller Performance Results 

The Lyapunov Hybrid Tracking Controller (LHPC), derived in Chapter 3, tied with 
the MHPC as the two best control algorithms for this phase of simulations. This controller 
uses an equilibrium glide phase for the first part of the trajectory and a Lyapunov Tracking 
exit phase using Lyapunov steepest descent techniques to steer the trajectory onto a target 
path. The target path selected is a 450 ft/sec constant altitude rate path which will lead the 
vehicle to the desired apocenter altitude.The transition velocity for switching from the 
equilibrium glide phase to the LHTC exit phase is varied using the technique of Chapter 2. 
The hybrid density estimation technique presented in Chapter 2 is used to define the de- 
sired path and to select the appropriate transition velocity. 

Testing this controller against the same perturbations considered earlier in this chap- 
ter produced excellent results. The results of injecting square wave density pulses into the 
nominal atmosphere, low density atmosphere, and high density atmosphere are summa- 
rized in Figs. 37 and 38, respectively. The results of varying L/D and entry flight path an- 
gle are illustrated in Fig. 39. 

This controller showed outstanding results for the density variations postulated for 
this set of simulations, showing practically no sensitivity to any of the perturbations. 
Again, however, the same caution presented in the MHPC performance results section 
should be repeated here: this simulated density is simply a function of altitude. The den- 
sity function for the outbound leg of the trajectory is identical to the density on the in- 
bound leg. The density estimation technique employed in this controller should excel 
under this condition. The robustness with respect to horizontal density variations must be 
evaluated to fairly generalize the evaluation of this (or any other) guidance scheme. 
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Lyapunov Tracking Controller Performance Results 

The LTC differs from the LHTC algorithm only in the density estimation technique. 
The LHTC employs the polynomial density estimator described in Chapter 2, while the 
LHTC uses the hybrid density estimation technique. The performance of this controller 
appears to be slightly degraded from the performance of the LHTC at about the same level 
that the performance of the MPC was worse than that of the MHPC. The performance is 
still acceptable and, as was stated in the analysis of the MPC’s performance, the strength 
of this density estimation technique is expected to surface when the outbound and inbound 
density functions differ. 

Figures 40, 41, and 42 illustrate the results of the square wave density pulses which 
perturb the nominal, low and high density atmosphere, respectively. Varying L/D and en- 
try flight path angle is depicted in Fig. 43. The worst performance noted during these sim- 
ulations using the LTC required 464 ft/sec to attain the desired orbit. This peak was 
caused by a +75% density pulse perturbing the high density atmosphere between 180,000 
and 200,000 ft. But again, even this worst case is considered to be acceptable. 

Selection of Controllers to Proceed 

The next stage of simulation was very intensive, requiring approximately sixty hours 
of computer time to fully test each controller. In an effort to limit this test matrix, only 
those controllers which showed promise of being able to handle the perturbations used in 
this chapter were to proceed to the next phase. The original plan was to select the two 
most promising controllers, and validate them. But, after analyzing the data presented in 
this chapter four controllers were selected for the next phase of testing. The four selected 
were the MHPC, the MPC, the LHTC and the LTC. The four selected are actually two 











94 


b) 


c) 

Fig. 43 Lyapi 
Ratio and Eni 



95 


control techniques, each employing two different methods of density estimation. All four 
of these controllers were able to handle the full range of testing performed during this 
phase without requiring more than 500 ft/sec to attain the desired orbit for any perturba- 
tion. The limitation of this set of simulations is that the inbound and outbound density 
functions are always the same. In the next chapter the performance of these four control 
algorithms will be investigated with different inbound and outbound density functions. 
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CHAPTER V 

DETERMINATION OF ROBUSTNESS LIMITS 

Determination of the robustness limits of the MHFC, MPC, LHTC and LTC, and se- 
lection of the most robust algorithm of these four is the goal of this chapter. Chapter 4 
showed that these four algorithms are all capable of handling extreme variations in the ve- 
hicle L/D, entry flight path angle and in the density function, provided the density function 
is a simple function of altitude. This chapter will examine the effect of more realistic den- 
sity functions which differ for the inbound and outbound legs of the trajectory. This will 
be accomplished by again injecting square wave density pulses into the density function, 
but this time the pulses will only perturb the outbound leg of the trajectory. In addition si- 
nusoidal variations in altitude and in vehicle range will be used to perturb the density func- 
tion. The control algorithms will also be tested using the actual density profiles measured 
by the Viking 1 and Viking 2 landers. 

To determine the robustness limits success and failure must first be defined. Because 
the vehicle has not been designed yet, the fuel budget for maneuvering the vehicle has not 
been defined. The definition of success and failure used here is somewhat arbitrary, 
though it is believed to be close to the actual definition. Success is defined as any aero- 
braking trajectory which requires 500 ft/sec or less of propulsive maneuvering (AV) to at- 
tain the desired final orbit. As in Chapter 4, AV is computed with three components, one 
propulsive maneuver applied at the atmospheric interface in the direction of the velocity 
vector to correct the apocenter altitude, a second at apocenter to raise pericenter and a 
third to correct any plane error. Because of the lack of a firm definition of vehicle charac- 
teristics a grey area has been defined. The grey area is any trajectory which requires be- 
tween 500 and 1,000 ft/sec. Any trajectory which requires between 500 and 1000 ft/sec to 
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attain the desired orbit will be referred to as a soft failure. Any vehicle design should cer- 
tainly carry enough fuel to perform 500 ft/sec of maneuvering to attain the desired orbit. 
Carrying additional fuel, however, to correct for a soft failure of the aerobraking control 
system reduces the advantage of aerobraking. The trajectories which terminate with soft 
failures would still be able to attain orbit, though not the desired orbit, using 500 ft/sec of 
propulsive maneuvering. This anomaly may result in some mission degradation, but not a 
complete mission failure. A hard failure is defined to be any trajectory which requires 
1,000 ft/sec or more of AVto attain the desired orbit. It includes any trajectories which fail 
to exit the atmosphere. By this definition, all four controllers considered in this chapter 
succeeded in all of the simulations performed in Chapter 4. 

Outbound Leg Square Wave Density Pulses 

The robustness test procedure begins by using square wave density pulses which per- 
turb the density function of the outbound leg only. The density during the descent into the 
atmosphere is either the nominal or a MARS-GRAM generated low or high density atmo- 
sphere model. After the altitude rate becomes positive, a square wave density pulse, simi- 
lar to those employed in Chapter 4 is used to perturb either a 10,000 or a 20,000 ft altitude 
band of the atmosphere. These pulses multiply the density predicted by the atmosphere 
model by 0.5, 0.75, 1.25, 1.5, 1.75 or 2.0 within the perturbed altitude band. The pulses 
are again referred to as -50%, -25%, +25%, +50%, +75% and +100% magnitude density 
pulses, respectively. As in Chapter 4, the pulses are moved in 10,000 ft altitude intervals, 
with the lower edge of the density pulse located between 100,000 and 290,000 ft. The per- 
formance is presented in Figs. 44 through 55 with AV plotted along the vertical axis while 
the magnitude of the pulse and the location of the lower edge are plotted on the other two 


axes. 
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MHPC Performance 

Performance of the MHPC in the nominal, low and high density atmospheres when 
the density of the outbound leg of the trajectory is perturbed by square wave density pulses 
is presented in Figs. 44, 45, and 46, respectively. In the first plot of each figure the pulse 
perturbs a 10,000 ft altitude band, while in the second plot the pulse perturbs a 20,000 ft 
altitude band. 

The MHPC produced many soft failures during this test sequence but no hard fail- 
ures were recorded. The 20,000 ft duration pulses produce worse performance than the 
10,000 ft duration pulses in almost all cases. The MHPC is very sensitive to large magni- 
tude (+50% and +75%) density pulses located below 180,000 ft in the nominal atmo- 
sphere, and 150,000 or 200,000 ft in the low or high density atmosphere, respectively. 
There is also a region of sensitivity caused by the -50% 20,000 ft density pulses. These 
latter regions are located at slightly higher altitudes and are not as severe as those caused 
by the larger magnitude pulses. In all of these plots there is a region at extremely low alti- 
tudes, where the pulses have little or no effect. This robust region occurs because these 
pulses are either located below the minimum altitude of the trajectory and the satellite nev- 
er flies in the perturbed atmosphere, or they are very near the minimum altitude of the tra- 
jectory and the satellite does not spend much time in the density fluctuations. 

There are two primary failure modes for these trajectories. When the large magni- 
tude density pulse affects the atmosphere in the altitude region where the satellite is in the 
equilibrium glide phase, the density filter is fooled into believing the entire atmosphere 
has higher density than that measured during the descent. The effect is to initiate the exit 
phase early, and predict a relatively high altitude rate for the exit phase. When the vehicle 
moves out of a high density region, there is a time lag before the density filter records the 
change. By the time the controller responds, the vehicle has moved even higher and the 
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vehicle has less control authority. As a result the vehicle leaves the atmosphere with too 
much energy, and overshoots the desired apocenter altitude. The second failure mode is 
caused by higher altitude density pulses. All of the failures caused by the -50% pulses ex- 
hibited this failure mode. The vehicle flies the equilibrium glide phase in the unperturbed 
atmosphere. After the vehicle initiates the exit phase it encounters the perturbed atmo- 
sphere. The large magnitude density pulses dissipate more energy than predicted, result- 
ing in a steeper exit phase, and in some cases, a lower apocenter than desired. Conversely, 
the small magnitude density pulses cause the vehicle to lose less energy than predicted and 
result in an apocenter altitude higher than desired. 

MPC Performance 

The MPC definitely performs better than the MHPC under these conditions. Again, 
the performance is presented in three figures, Figs. 47, 48, and 49, with the first figure 
showing results from the nominal atmosphere, the second from the low density atmo- 
sphere and the third from the high density atmosphere. 

The 10,000 ft density pulses has almost no effect on the performance of this control 
algorithm. Even the 20,000 ft pulses produce reasonably good results. There is only one 
soft failure noted during this test sequence and two very near failures for the MPC. All 
three of these events are caused by 20,000 ft +100% density pulses perturbing the low den- 
sity atmosphere. The pulse between 120,000 and 140,000 ft requires a AV of 584 ft/sec, 
while the pulse 10,000 ft higher (between 130,000 and 150,000 ft) requires 498ft/sec. The 
pulse between 140,000 and 160,000 ft requires just 458 ft/sec, but the one between 
150,000 and 170,000 ft requires 492 ft/sec. These high AVs are caused by the same failure 
modes as described above with the pulses having lower edges at 120,000 and 130,000 ft 
causing the density estimator to overreact and force the vehicle to exit with too much ener- 
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gy, while the pulse at 150,000 ft causes the vehicle to lose more energy than planned. 
Overall, the polynomial density estimator used in the MFC shows improvement over the 
hybrid density estimator of the MHPC. As suggested in Chapter 4, the promise of this 
density estimation technique is realized when the inbound and outbound density functions 
are different. 

LHTC Performance 

Performance of the LHTC is mixed. For the majority of density perturbations calcu- 
lated here this controller performs better than either the MHPC or the MPC. Yet there are 
a few isolated instances where the controller performs extremely poorly. The controller 
even produces two hard failures. Presentation of this controller’s performance follows the 
same format as before with the nominal atmosphere results in Fig. 50, the low density at- 
mosphere results in Fig. 51 and the high density atmosphere results in Fig. 52. 

One noteworthy aspect of the Lyapunov Tracking exit phase is that it almost always 
commands either full lift up or full lift down. When the vehicle is exactly on the desired 
path, the commanded bank angle chatters between ±15° and ±165° (commanded bank an- 
gles less than 15 or greater than 165° are allowed only when the orbit plane error is less 
than .03°). Of course, the vehicle roll rate and roll acceleration limits prevent the vehicle 
from oscillating too wildly. But, when the vehicle is not on the desired path the control 
system will command near full lift up, or full lift down to approach the trajectory. This 
feature allows the vehicle to respond more quickly than it does for the predictor corrector 
algorithms to pull the vehicle back onto the desired path. However, when the desired path 
is computed poorly because of a poor density estimate, the control system still responds by 
commanding full lift to approach the computed path as rapidly as possible. This controller 
also suffers from the same problems with the density estimator as the MHPC. This phe- 
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nomena causes the two hard failures seen in Fig. 50b and Fig. 52b. The onboard acceler- 
ometer measurements records high drag measurements while the vehicle decelerates in the 
region of high density caused by the density pulse. The density filter predicts that, be- 
cause of the measured high drags, the remainder of the atmosphere is also of higher densi- 
ty than that recorded during the descent. The effect is to predict higher energy loss due to 
this drag than actually occurs. This inaccurate prediction causes the control system to ini- 
tiate the exit phase earlier than desired and to command a path which climbs out of the at- 
mosphere at relatively high speed. As noted above, the Lyapunov optimal control solution 
pulls onto this path as rapidly as possible. By the time the satellite moves out of the high 
density region, and the density filter recognizes the change, the satellite is too high and is 
at a velocity which is too high too allow recovery. The result wis a post aerobraking apo- 
center altitude of 836 nm for the hard failure in Fig. 50b and 1,165 nm in Fig. 52b. Again, 
the target apocenter is 270 nm. 

The Lyapunov Tracking exit phase, in spite of the two hard failures described in the 
previous paragraph, seems better able to cope with these density estimation problems than 
the predictor corrector algorithms. The rapid response of the vehicle is largely advanta- 
geous. The nature of the Lyapunov control system, though its agressive technique causes 
the two hard failures discussed above, is still desirable. As accelerometer measurements 
are taken and the density filter is continually updated, the desired path varies. The rapid 
response of the Lyapunov control law helps track this moving path as long as the vehicle 
has enough aerodynamic control authority to respond. In the LTC controller simulation 
results which follow, the effect of combining the polynomial density estimation technique 
with the fast response of the Lyapunov control scheme is further explored. 
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LTC Performance 

The LTC performance is presented in Figs. 53, 54, and 55. The 10,000 ft density 
pulses have almost negligible effect on the performance of this control algorithm, as is the 
case with the MPC. Also, the 20,000 ft pulses produce reasonably good results with the 
LTC. There are only three soft failures observed during this test sequence for the LTC. 
All three of the failures are caused by the same 20,000 ft +100% density pulses perturbing 
the low density atmosphere which caused the soft failure and the two other near failures in 
the MPC performance. The pulse between 120,000 and 140,000 ft requires a AV of 542 ft/ 
sec, while the pulse 10,000 ft higher (between 130,000 and 150,000 ft) requires 505 ft/sec. 
The pulse between 140,000 and 160,000 ft does not result in a soft failure, requiring 493 
ft/sec, but the one between 150,000 and 170,000 ft does, requiring 533 ft/sec. These fail- 
ures are caused by the same failure modes as described earlier. The pulses located at 
120,000 and 130,000 ft cause the density estimator to overreact, forcing the vehicle to exit 
with too much energy, while the pulse at 150,000 ft causes the vehicle to lose more energy 

than planned. 

Overall, though, the polynomial density estimator used in combination with the Ly- 
apunov control scheme shows excellent performance. The performance of this control al- 
gorithm during this evaluation sequence very nearly parallels that of the MPC. 


Sinusoidal Density Variations 


The next set of simulations involves perturbing the density function with sine waves. 
Sine waves in altitude and sine waves in range are used. These sine waves vary in ampli- 
tude (K a ), wavelength ( X ) and phase angle (4>). The sine wave perturbations in altitude 


took the form 
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while those in range took the form 
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The range of amplitudes used included 0.1, 0.25 and 0.5 for both forms of penurba- 

tions. The wavelengths selected for the altitude variations included 1000, 2000, 5000, 

10,000, 20,000, 50,000, 100,000, 200,000 and 500.000 feet. For the variations in range 

the wavelengths selected included 1, 5, 10, 20, 50, 100, 200, 500. 1000, 2000 and 5000 

In . n , 

nautical miles. In both cases, the phase angles included zero through — in - increments. 


The sinusoidal variations with amplitude of 0.1 appears to be very much in line with 
the actual density profiles measured by Viking 1 and Viking 2 landers during their descent 
through the Martian atmosphere (Figs. 7 and 8). All of the results with an amplitude of 0. 1 
for both forms of the sinusoidal variations and all four controllers examined in this chapter 
are successful. The highest AV requires 490 ft/sec, but the vast majority of the trajectories 
(over 99%) require less than 400 ft/sec. Only 13 of the 1920 trajectories tested with a 0.1 
amplitude sine wave density variation need more than 400 ft/sec of AV Likewise, the re- 
sults generated using 25% and 50% amplitude sine waves in altitude are almost as benign 
as the 10% results. All of the trajectories which used 25% amplitude sine waves in alti- 
tude are successful. Of the 864 trajectories checked using 50% amplitude sine waves in 
altitude, none result in hard failures and only 8 produce soft failures. Of these, only 3 de- 
mand more than 600 ft/sec, with the worst calling for 732 ft/sec. A complete breakdown 
of these failures is presented in Table 2. Because these results generated using 10% sine 
waves in altitude and range and 25% sine waves in altitude are all successful, and the 8 
soft failures generated using 50% sine waves in altitude are adequately described in Table 
2, they are notpresented graphically. It is interesting to note that 3 of the 4 soft failures 
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which occurred with controllers using the hybrid density estimation technique, have wave- 
lengths of 1000 ft, while the fourth has a wavelength of 5000 ft. The 1000 ft wavelength 
sine waves in altitude seem to be corrupting the stored density data used in the density es- 
timation process. Thess data are stored at 1000 ft altitude intervals. Though the density 
filter should be able to compensate for this eventuality, it does not appear to do so well 
enough to prevent these failures. Three of the four failures which occurred with control- 
lers employing the polynomial density estimation technique have wavelengths of 20,000 
and 50,000 ft. Shorter wavelengths tend to have a cancelling effect, with the additional 
drag of high density regions being offset by the lower drag of low density regions. Longer 
wavelengths are easy for the sixth order polynomial to follow, provided there are no more 
than five extremes in the density function. The problem with the 20,000 and 50,000 ft 
wavelength sine waves is they do not oscillate fast enough to cancel high density regions 
against low density regions, yet they still have six to fifteen complete sine waves, with 
twelve to thirty density extremes in the aerobraking region. This variation is simply more 
than a sixth order polynomial can follow. The final failure is caused by an excessive orbit 
plane error (wedge angle) at exit. 

25% and 50% Sine Waves in Range 

The 25% amplitude sine wave density perturbations, which use vehicle range from 
entry as the argument to the sine function, are probably the truest measure used here to test 
controller robustness in the presence of a realistic worst case Martian atmosphere. These 
perturbations are of somewhat higher magnitude than the perturbations measured by the 
Viking 1 and Viking 2 landers, but, most probably, the Viking 1 and Viking 2 landers did 
not sample the worst case atmospheric perturbations. Though the amplitude of the per- 
turbing sine wave is increased to 50% for the simulations, the probability is extremely low 
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that the Martian atmosphere ever experiences high frequency oscillations in density with- 
amplitudes this large. 

The performance of the MHPC when tested against the 25% amplitude perturbations 
is presented in Figs. 56a, 57a, and 58a for the nominal, low density and high density atmo- 
spheres respectively. Figures 56b, 57b, and 58b present similar results when the amplitude 
of the perturbing sine wave increases to 50%. Similarly, Figs. 60 and 61 present the re- 
sults for the MPC, while the results for the LHTC are in Figs. 62, 63, and 64 and the LTC 
results are depicted in Figs. 65, 66, and 67. 

The 25% amplitude perturbations are significant enough to cause problems in some 
of the trajectories. Though they do not induce any hard failures, there are many soft fail- 
ures. The 50% amplitude perturbations are severe enough to cause several hard failures 
for all of the controllers except the LTC. The 25% and the 50% amplitude sine waves are 
each used to simulate 264 perturbed atmospheres for each controller (11 wavelengths x 8 
phase angles x 3 base atmospheres). Of these 264 trajectories tested using the MHPC and 
the 25% amplitude sinusoidal variation, six trajectories result in soft failures. Six trajec- 
tories also resulted in soft failures when the MPC controller is used, though they are not 
the same six perturbations. The LHTC has four soft failures while the LTC only has two. 
When the amplitude of the perturbing sine wave is increased to 50%, the MHPC had four- 
teen hard failures, the MPC had ten and the LHTC had fourteen. These three controllers 
also experienced many soft failures during these simulations. The LTC did not result in 
any hard failures, but it did produce twenty-nine soft failures. 

All of these failures are the result of exit phase failures, which are in turn attributable 
to related directly to density estimation difficulties. Nonetheless, the equilibrium glide 
phase is robust enough to keep the vehicle in the atmosphere and prevent a skip out for all 
of these trajectories. None of the trajectories fail to exit the atmosphere, although some of 
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them barely do so. The problem with all of the failures centered around the inability of the 
density estimation technique to adequately predict the density function and thus the 
amount of drag expected by the controller for the duration of the trajectory. Both density 
estimation techniques appropriately ignored the high frequency density variations (those 
with wavelength less than 10 nm). These oscillations occur so quickly that the high and 
low density regions have a cancelling effect. 

The hybrid density estimator shows increased sensitivity to wavelengths of 20 to 200 
nm. The polynomial density estimator, on the other hand handles these wavelengths very 
well. It is the 500 to 2000 nm wavelengths which produce problems for this estimator. 
These sensitivities to different wavelengths are easy to understand. The hybrid density es- 
timation technique uses the density filter to adjust its estimate for the entire atmosphere 
based on current density measurements. The long wavelength sine waves have the same 
effect as a slowly increasing or decreasing density bias during the trajectory. The density 
filter of the hybrid density estimator is able to sense this slow drift and appropriately adjust 
the measurements taken during descent to compensate for the drift. The wavelengths 
which give the hybrid density estimator trouble are those which perturb a portion of the at- 
mosphere and then reverse that perturbation fast enough to confuse the density filter but 
not fast enough to have a cancelling effect. The polynomial density estimator, on the other 
hand, fits the sixth order polynomial in altitude to the normalized density function. This 
density estimation technique remembers the density which was measured at the various al- 
titude intervals. It takes the most recent density measurement and adds this information to 
the knowledge base and fits a smooth polynomial curve through the data. When the local 
density is biased, but then that bias reverses later in the trajectory, as it does when the in- 
termediate wavelength sine waves perturb the atmosphere, this density estimation tech- 
nique excels. But, when the density function is monotonically increasing or decreasing 
during the trajectory, as is the case for the longer wavelength sine waves, this estimation 
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technique does not respond fast enough. An attempt to place more weight on the most re- 
cent data should help, but attempts to do so made the oldest data obsolete; that is, the high- 
er altitude densities, with the density estimator sometimes missing the density at exit by an 
order of magnitude or more. Clearly, the development of better density estimation mod- 
ules deserves further study. 

Overall, the Lyapunov control scheme performed better than the predictor corrector. 
The rapid response of the Lyapunov tracking exit phase compensates well for slowly de- 
veloping density estimates. The polynomial density estimator also performs better than 
the hybrid density estimator, as it most clearly seen from the LTC results. The LTC kept 
AV below 500 ft/sec for all but two of the 25% amplitude sine wave perturbed atmo- 
spheres, and those two only required a AV of 509 or 577 ft/sec. Additionally, the 50% 
amplitude trajectories are all completed with AV below 1000 ft/sec. The LTC also copes 
with the square wave density pulses, both those presented in Chapter 4 which perturbed 
the entire atmosphere and those of this chapter which only effect the outbound leg of the 
trajectory. Since the LTC required less than 500 ft/sec for all of the trajectories tested in 
Chapter 4, and responded better than any of the other controllers to all the robustness tests 
of this chapter, the LTC is the most robust aerobraking controller examined. 
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CHAPTER VI 

CONCLUSIONS AND RECOMMENDATIONS 


Conclusions 


The Analytic Predictor Corrector algorithm selected as the control algorithm for the 
AFE is generally a robust control algorithm, especially with respect to large scale density 
variations. The algorithm is fairly robust to short period density variations, but does dem- 
onstrate a definite sensitivity to variations in the entry flight path angle and vehicle lift and 
drag coefficients. These sensitivities are due, in large part, to the fixed transition velocity 
employed to switch the control algorithm from the entry phase to the exit phase and the 
rather simplistic density estimation scheme used. It is necessary to increase the K- term 
in the equilibrium glide phase to prevent rapid large scale density variations from causing 
a premature exit from the Martian atmosphere. 

The Energy Controller is slightly more robust than the APC to variations in the entry 
flight path angle and to uncertainty in vehicle lift and drag coefficients. It is also robust to 
large scale density variations. However, short period density variations were unacceptable 
to this control algorithm and the increased trajectory loads caused by the EC led to its ear- 
ly dismissal from the list of potential control algorithms. 

The Numerical Gradient technique, and then the Conjugate Gradient technique are 
used to compute idealized optimal (minimum AVI trajectories. It was hoped that these 
methods could be adapted as an on-board control algorithm. But. these algorithms require 
about two orders of magnitude more computational time than the APC or EC to generate a 
solution. Additionally, the optimization technique assumes all pertinent density and vehi- 
cle lift and drag characteristics are known precisely. The trajectories produced by these 
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optimization techniques fly the exit phase using the full lift available to remain in the at- 
mosphere. Any decrease in density from that modeled in the optimization process allows 
the vehicle to exit early with too much velocity. These algorithms, with the current perfor- 
mance index, do not produce robust trajectories even if they are able to compute a solution 
fast enough to control the satellite in real time. A more general performance index which 
seeks to minimize AV while retaining robustness and also reducing control activity should 
be sought if these techniques are to become practical. 

The modifications proposed to the APC to produce the MHPC and MPC convert that 
algorithm into a robust control algorithm capable of guiding the aerobraking trajectory to 
near minimum AV exit state for most of the perturbations considered 

s mentioned before, it is necessary to increase K- for the equilibrium glide phase to 
prevent a premature exit from the atmosphere. But in addition, the change to the more 
computationally straight forward and efficient exit phase, combined with the better density 
estimation techniques and the variable transition velocity, made significant improvements 
to overall robustness of the control algorithms. Between the MHPC and the MPC, the 
MPC responded better on the whole to the perturbations examined here. There were two 
areas where the MHPC did slightly better than the MPC. The first situations occurs when 
density is simply a function of altitude and the entry and exit density functions are identi- 
cal. The probability of such a coincidence is rather low, but the MPC is still able to handle 
these situations well (though not as well as the MHPC), without producing any failures. 
The second area is when the large amplitude sinusoidal variations, which used range from 
entry as the argument to the sine function, had wavelengths between 500 and 2000 nm. 
This possibility is still a concern and leads to several of the recommendations that follow. 

Overall, however, the MPC reacted more appropriately to realistic perturbations than did 
the MHPC. 
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The Lyapunov Steepest Descent Control algorithm is implemented, but its inability 
to compensate for varied energy depletion rates due to density variations, variations in en- 
try flight path angle or vehicle drag coefficient made this algorithm unusable. However, 
when the Lyapunov control algorithm is recast as a tracking controller designed to follow 
a reference trajectory, it shows much more promise. The algorithm still has trouble exiting 
with just the right amount of energy to achieve the desired apocenter altitude and produces 
peak trajectory loads higher than those of the predictor corrector algorithms. To cure the 
first ailment a scheme to vary the gain values in the Lyapunov function is developed, while 
the second is fixed by employing the equilibrium glide entry phase and using the Ly- 
apunov Tracking Controller as an exit phase. 

With the two density estimation techniques developed for the MHPC and the MPC 
used to define the reference trajectory, and the transition velocity from entry to exit phase 
computed as for the predictor correctors, the LHTC and LTC performed extremely well. 
The performance of the LHTC and LTC essentially mirror that of the MHPC and MPC, re- 
spectively. Generally, the strengths of the MHPC are the strong points of the LHTC, while 
they share common weaknesses as well. Perturbations which cause problems for the MPC 
are also likely to cause problems for the LTC. In most cases, the problems are initiated be- 
cause the density estimation technique is unable to follow a specific perturbation. The Ly- 
apunov tracking algorithm, with its more rapid response, compensates better and produces 
exit states which require less AV than the predictor correctors. There were a few notable 
exceptions where the rapid response moved the vehicle into a less dense region too rapidly 
resulting in loss of control authority and an exit state with too much energy. But, predom- 
inantly, the Lyapunov trackers performs better than the predictor correctors. As in the pre- 
dictor corrector analysis, the polynomial density estimation technique works better than 
the hybrid density estimation technique. Overall, the LTC performs better than the 
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LHTC, MPC or MHPC and is the recommended control algorithm for performing an in- 
terplanetary aerobraking maneuver at Mars. 

Recommendations 

Based on these conclusions, the following recommendations are made: 

1) Robustness to density variations should be a prime issue in selecting the control al- 
gorithm for the aerobraking phase of the MRSR. This characteristic must be consid- 
ered along with decisions such as entry velocity, vehicle lift requirements, ballistic 
coefficient, or navigational accuracy requirements. 

2 j The expected wavelengths and maximum amplitude of the short period density oscil- 
lations in the Martian atmosphere should be characterized. The nature of these short 
period oscillations should be determined. It would be beneficial in designing a den- 
sity estimation technique to know if the short period density wave structure is prima- 
rily horizontal or vertical in nature, or a predominantly time-varying function. 

3) Once the frequency of the expected density variations is determined, the density esti- 
mation technique employed in the aerobraking control system should be tuned to re- 
spond to the most likely frequencies which may perturb the trajectory, while 
ignoring those which have minimal effect on the trajectory. 

4) A higher order density estimator, perhaps using Tschebechev polynomials or Leg- 
endre polynomials to bypass the numerical difficulties of a higher order polynomial 
in altitude, should be examined. It may also be desirable to fit a second function, in 
terms of arc length, or time, or range to the density function, especially if a monoton- 
ically increasing or decreasing density function is predicted. 
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5) The density estimation technique should be adjusted to use all available knowledge 
of the Martian atmosphere, including any knowledge of dust storms, solar flares, or 
solar heating of the atmosphere along the intended trajectory. 

6) The LTC should be tested using higher entry velocities and different vehicle lift and 
drag characteristics or ballistic coefficient, as well as different target orbits to deter- 
mine its suitability for controlling some of the other mission scenarios proposed for 
MRSR. The fast trip manned precursor mission is clearly a candidate. Also, trading 
off nominal performance for robustness by varying the exit phase altitude rate should 
be studied. 

7) A statistical method of evaluating controller performance should be developed based 
on the probability of various atmospheric perturbations occurring. This method may 
extend further to include the probability of variations in entry conditions or vehicle 
aerodynamic characteristics. 

8) A new performance index should be developed which will minimize AV while re- 
taining a level of robustness. With this new performance index, the calculus of vari- 
ations optimization techniques should be revisited in an attempt to construct a 
controller which computes a truly optimal solution. 
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APPENDIX A 

IDEALIZED MINIMUM AV OPTIMAL SOLUTION 

A numerical gradient technique was employed to determine the minimum AV solu- 
tion for a nominal Martian aerobraking maneuver 40 . The MRSR mission scenario calls 
for the aerobraking maneuver to reduce the vehicle velocity relative to the planet using 
aerodynamic drag and then to exit the atmosphere on an elliptical intermediate orbit. A 
series of propulsive maneuvers are then performed to transfer the vehicle from the inter- 
mediate orbit to the desired final orbit. The total AV required to transition from the inter- 
mediate orbit to the desired orbit is determined by the vehicle’s atmosphenc exit velocity 
vector and is a good measure of control system performance. The open loop solution pre- 
sented here assumes that initial conditions as well as all pertinent vehicle and atmospheric 
properties are known precisely. Limits are not placed on trajectory loads. Robustness to 
atmospheric dispersions is not considered in computing this optimal solution. This solu- 
tion produces the minimum AV attainable to transition from the post aerobraking interme- 
diate orbit to the desired final orbit for a given atmosphere, vehicle and entry condition and 
is used as a benchmark to evaluate the performance of the feedback controllers. 


Equations of Motion 

The formulation begins with the equations of motion. The equations of motion were 
presented in Chapter III but are repeated again here for completeness. 


dr 

dt 


Q = Vsiny 
dt 


(86) 
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dt 
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( 88 ) 


Equation (86) is simply the radial velocity in terms of the inertial velocity and flight 
path angle. Equation (87) gives the time rate of change of velocity in two parts: 1 ) the ve- 
locity loss rate due to aerodynamic drag and 2) the change in velocity due to gravitational 
acceleration (the inertial component). Similarly, Eq. (88) is the time rate of change in the 
flight path angle also composed of two parts: 1) the change in flight path angle due to the 
component of aerodynamic lift in the vertical plane and 2) the change in flight path angle 
due to gravitational acceleration (the inertial component). The control variable <t>, the 
bank angle, determines the amount of lift exerted in the vertical plane to bend the trajecto- 
ry and change the flight path angle. 


Nondimensional State Variables 

Dimensionless state variables are introduced: 



X 1 


h/h 
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X = 

X 2 
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V/iJp/R) 
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along with a dimensionless time variable x 

x = ( t/h e )JpVR . 


(89) 


(90) 


The equations of motion may now be written: 
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where a = p/p 0 = exp [ (— (A — h 0 ) ) /hS ] , A - (p Q Sh e Cj) / (2m) 
B = (p 0 S/i e C D )/(2m)andc = 


The Performance Index 


To minimize total AV required to transition to the desired orbit it is sufficient to min- 
imize the exit flight path angle provided the apocenter of the intermediate orbit equals the 
desired apocenter. This procedure maximizes the pericenter of the post-aero orbit. Two 
terminal constraints are employed. The first requires the final altitude to be the atmospher- 
ic interface altitude and the second fixes the intermediate orbit apocenter. The cost func- 
tion is therefore the exit flight path angle (7 = y x ) and the goal is to minimize the cost 
function subject to 
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( 95 ) 
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To derive Eq. (95), set the orbital angular momentum at exit equal to the angular mo- 
mentum at apocenter. 


h = r x V x cosy x = (96) 

From this equation solve for the velocity at apocenter in terms of the terminal radius, 
velocity, flight path angle and the radius of apocenter. 


V « = 


r X V X C0S Yr 


(97) 


Equate the orbital energy at exit to that at apocenter, using the expression for veloci- 
ty at apocenter from above 


(f JC V x cosy v ) 2 


2ri 
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r 



(98) 


Obtain Eq. (95) after some algebra and after replacing the physical state variables 
with the nondimensional variables given in Eqs. (89) and (90). 


The Numerical Gradient Technique 

This problem is solved using a first order numerical gradient procedure. To formu- 
late the optimal control problem begin with the performance index. In general terms this 
index may be written 


J = t(x f ) +j'f o {Hx*,u*)}dt 


(99) 
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The performance index is then augmented with penalty functions to impose the ter- 
minal constraints and the equations of motion. 


J a = $( Xf ) + [v] T {V(x f )} + {L (x * , u * ) + [X] T \f(x*,u*) -x]}dt (100) 

For this problem [v] is a 2 x 1 column matrix of constants, (y (*/) } is a 2 x 1 
column matrix of terminal constraints given by Eqs. (94) and (95) and [X] is a 3 x 1 time 
varying matrix of Lagrange multipliers or influence functions, {x} = {f (x, u) } are the 
3 first order differential equations of motion Eqs. (91), (92) and (93). { u} is the control 
variable O. To customize this general augmented cost function for the problem at hand 
delete L (x * , u ) since there is not an integral term in our performance measure, and sub- 
stitute jj- for <{) (Xf) to obtain 

J a = y f + [v] 7 ’(vU/)> +I^W r l f(x*,u)-x]}dt (101) 

The numerical solution process begins with a guess of the control time history. The 
values for the state variables are computed from initial conditions and then integrated for- 
ward in time using this postulated control time history. Differential equations for the 
Lagrange multipliers, which are developed later, are used to integrate [X] backward in 
time beginning with the value of [X.] computed at tj. A new control time history is de- 
rived by setting the first variation in the augmented cost function to zero. The process is 
repeated until the terminal constraints are satisfied to within an acceptable tolerance. 

The first variation in J a is formed 


bJ a = ^8^+ [v] r ^8x /+ {'( W r ^) 5* + < W T |;> s “ - W T&i > d ‘ ■ < 102 > 

/ f 
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Integrate [X] T bx by parts to obtain 

-ft ( [X] T bx) dt = -[X f ] T bx f + J *{ ( [X] T bx) dt (103) 


Substitute Eq. (103) into Eq. (102) to obtain 
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(104) 


Define a new Lagrange multiplier 


[\] T = [X J ] T + [v ] T [X 1 ] T 


(105) 


[X^] is the 3x1 column of Lagrange multipliers normally used to impose the 
equations of motion while [X 1 ] is dimensioned 3x2 and contains additions to the 
Lagrange multipliers which arise from the terminal constraints. The first variation in J Q 
may now be written 
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(106) 
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The performance index is minimized by setting the first variation in J a equal to zero. 
Choose differential equations for the Lagrange multipliers so that the coefficient of bx 
goes to zero to obtain the two equations 

[i'] r = -[V'] 7 '| (107) 


and 
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Integrate the costates backwards in time using Eqs. (107) and (108) with the bound- 
ary conditions obtained by requiring the coefficient of 8.xy in Eq. (106) to be zero 
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Since the coefficients of 8 jc and Sxj have been set to zero, the first variation of J a re- 


duces to 
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Defining two new variables 
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is a scalar, while A 11( is a 2 x 1 column matrix and 
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= jo 0 -Aax 2 sinOj . 
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To make J Q as small as possible, choose the variation in the control 8m to be 

8m = -K [A^ + [v] t A X¥ ] T (127) 

AT is a scalar weight which fixes the relative importance placed on minimizing the 
cost function versus satisfying the terminal constraints. A value of 200 for K places suffi- 
cient weight on the cost function and still allows the terminal constraints to be satisfied 
within an acceptable tolerance. By substituting Eqs. (127) and (125) into Eq. (123) obtain 


{8\j/} = [A^+ [v] ^A^] T dt. (128) 

Again, introduce two additional variables 

{g} [\J [A*] 7 ’** (129) 

[Q] s j! 7 [A v ] [\] T dt (130) 

Substitute these variables into Eq. (128) to obtain 

{8\j/} = -K [g + G[v]] . (131) 


To drive {8y} to zero, choose {8y} = — { \py} , where {vj/y} is the value of the 
terminal constraints, computed after integrating the state equations forward, solving Eq. 
(131) for {v} 


{V} = 


(132) 
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Use this value for {v} in Eq. (127) to obtain (5 «} . The control update is then 
computed 


ne = t u old > + • 


(133) 


Since the final time is free, we must also minimize 


"“Sr,= [*+^8, 

dtf f [dtf 


37 / ' dt f \°T 


(134) 


Now let 


5 '/ = -i 


3(1) 7’3\it' 

4. vr i — 

L */ 


(135) 


Replace 4> with y f and y with the expressions given in Eqs. (94) and (95) to obtain 



(136) 


Use the new control time history {Eq. (133)}, along with the change in t f {Eq. 
(136)}, to again integrate the state equations forward. Compute the terminal value of the 
Lagrange multipliers and integrate these backwards in time, then recompute the control 
time history. This process is repeated until the terminal constraints are satisfied within an 
acceptable error bound. The final apocenter altitude is required to be within 5 nm of the 
target value while the terminal altitude is required to be within 25,000 ft of the defined at- 
mospheric interface altitude. Apocenter errors of 5 nm require very little AV to correct 
and are attainable using this optimization method although thirty or more iterations may 
be required to converge this closely. The 25,000 ft terminal altitude error band was chosen 
because the aerodynamic effects decrease exponentially with altitude and are almost negli- 
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gible at altitudes above 300,000 ft. To converge closer than 25,000 ft to the selected atmo- 
spheric interface altitude of 410,105 ft (125 km) requires many more iterations. 
Furthermore, the terminal altitude generally converges toward the target altitude from 
above. 


Conjugate Gradient Projection Method 

The conjugate gradient projection method was employed to speed convergence of 
this problem 40, 41 . The gradient obtained in Eq. (127) was again used in this method to 
compute the search direction for correcting the control variable. However, after the first 
control update the previous search direction is used in conjunction with the computed gra- 
dient to give the problem near second-order convergence characteristics. The procedure 
follows. First, compute the gradient direction using Eq. (127) 

8i = -8 u (137) 

where the i subscript denotes the z'th iteration of control updates. Next, compute the search 
direction 


s = — P. 
i 6 i 


for the first iteration, while for subsequent iterations 


(138) 


s i = -£; + 


( £,•> 8i ) 


S : 




(139) 


where (a, b ) is the inner product of a and b. 

Once the search direction is determined, it is necessary to properly scale the magni- 


tude of the correction. 





